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ABSTRACT 


Swarm  attacks  are  of  great  concern  to  the  U.S.  Navy  as  well  as  to  navies  around 
the  world  and  commercial  ships  transiting  through  waters  with  high  volume  of 
marine  traffic.  A  large  group  of  hostile  ships  can  hide  themselves  among  various 
other  small  ships,  like  pleasure  crafts,  fishing  boats  and  transport  vessels,  and 
can  make  a  coordinated  attack  against  a  High  Value  Unit  (HVU)  while  it  passes 
by.  The  HVU  can  easily  be  overwhelmed  by  the  numbers  and  sustain  heavy 
damage  or  risk  being  taken  over. 

The  objective  of  this  thesis  is  to  develop  heuristic  algorithms  that  multiple 
defenders  can  use  to  intercept  and  stop  the  advances  of  multiple  attackers.  The 
attackers  are  in  much  larger  numbers  compared  to  the  defenders,  and  are 
moving  in  on  a  slow  moving  HVU.  Pursuit  guidance  laws  and  proportional 
navigation  (PN)  guidance  laws,  commonly  used  in  missile  guidance  strategies, 
are  modified  to  be  used  by  the  defenders  to  try  intercepting  attackers  that 
outnumber  them. 

Another  objective  is  to  evaluate  the  effectiveness  of  the  heuristic 
algorithms  in  defending  the  HVU  against  the  swarm  attack.  The  probability  that 
the  HVU  survives  the  swarm  attack  will  be  used  as  a  measure  of  effectiveness  of 
the  algorithms.  The  impact  of  various  parameters,  like  the  number  of  defenders 
and  the  speed  of  defenders,  on  the  effectiveness  of  the  algorithms  are  also 
evaluated. 
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I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

Since  the  attack  on  U.S.S.  Cole  in  2000  and  the  U.S. -Iranian  naval  dispute 
in  the  Strait  of  Hormuz  in  2008,  awareness  for  asymmetric  warfare  has 
increased.  Swarm  attacks,  in  particular,  are  of  great  concern  to  the  U.S.  Navy  as 
well  as  to  navies  around  the  world  and  commercial  ships  transiting  through 
waters  with  high  volume  of  marine  traffic.  A  large  group  of  hostile  ships  can  hide 
themselves  among  various  other  small  ships,  like  pleasure  crafts,  fishing  boats 
and  transport  vessels,  and  can  make  a  coordinated  attack  against  a  High  Value 
Unit  (HVU)  while  it  passes  by.  The  HVU  can  easily  be  overwhelmed  by  the 
numbers  and  sustain  heavy  damage  or  risk  being  taken  over. 

The  NPS  thesis  by  Tiwari,  2008  [1]  presented  the  capability  gap  in 
defending  against  such  a  swarm  attack.  This  thesis  aims  to  study  defensive 
strategies  that  can  be  used  to  defend  against  swarm  attacks,  and  evaluates  the 
mission  effectiveness  of  these  strategies  using  a  cost  function  that  can 
realistically  represent  a  multi-agent  engagement  scenario. 

B.  OBJECTIVES 

The  objective  of  this  thesis  is  to  develop  heuristic  algorithms  that  multiple 
defenders  can  use  to  intercept  and  stop  the  advances  of  multiple  attackers.  The 
attackers  are  in  much  larger  numbers  compared  to  the  defenders,  and  are 
moving  in  on  a  slow  moving  HVU.  Pursuit  guidance  laws  and  proportional 
navigation  (PN)  guidance  laws,  commonly  used  in  missile  applications,  are 
modified  to  be  used  by  the  defenders  to  try  intercepting  attackers  that 
outnumbers  them. 

Another  objective  is  to  evaluate  the  effectiveness  of  the  heuristic 
algorithms  in  defending  the  HVU  against  the  swarm  attack.  The  probability  that 
the  HVU  survives  the  swarm  attack  will  be  used  as  a  measure  of  effectiveness  of 
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the  algorithms.  The  impact  of  various  parameters,  like  the  number  of  defenders 
and  the  speed  of  defenders,  on  the  effectiveness  of  the  algorithms  are  also 
evaluated. 

C.  LITERATURE  REVIEW 

Various  existing  literature  related  to  multi-agent  engagement  scenarios 
usually  focus  on  search  algorithms.  Chung  et  al,  2011  [2]  looks  at  optimal 
detection  of  underwater  intruder  in  a  channel  using  Unmanned  Underwater 
Vehicles.  Royset  &  Sato,  2010  [3]  considers  route  optimization  for  multiple 
searchers  to  look  for  one  or  more  probabilistically  moving  target(s).  In  these  two 
papers,  there  is  only  one  target  with  multiple  defenders  searching  for  it.  Jang  & 
Tomlin,  2005  [4],  Shaferman  &  Oshman,  2009  [5],  Shin,  2011  [6]  investigates 
various  cooperative  guidance  strategies  in  multi-agent  engagement  scenarios, 
but  are  mostly  limited  to  missile  guidance  and  considers  scenarios  with  only  3  or 
4  targets,  defended  with  an  equal  number  of  defenders.  Rozen,  2009  [7],  studies 
the  detection,  recognition  and  interception  of  multiple  targets  using  an  interdiction 
force  of  a  UAV  and  a  navy  vessel,  but  the  number  of  targets  are  small  and 
moves  at  random  and  cannot  be  considered  as  a  swarm.  In  addition,  mission 
effectiveness  is  based  only  on  the  number  of  targets  considered  intercepted 
under  conditions  that  are  not  reflective  of  real  life  situations. 

This  thesis  explores  strategies  that  naval  ships  can  use  to  engage 
attackers  that  greatly  outnumbers  themselves  effectively,  and  uses  a  cost 
function  that  considers  the  hit  rate  of  the  defenders  against  the  attackers  and  the 
attackers  against  the  HVU,  to  work  out  the  mission  effectiveness  by  evaluating 
the  survival  rate  of  the  HVU. 

D.  THESIS  ORGANIZATION 

In  the  next  chapter,  the  setup  of  the  simulation  is  explained,  describing  the 

parameters  used  to  simulate  the  motion  of  the  HVU,  the  attackers  and  the 

defenders.  The  derivation  of  the  cost  function  used  to  evaluate  the  mission 

effectiveness  of  the  algorithm  will  also  be  shown.  A  test  matrix  is  formulated  to 

2 


investigate  the  effects  of  adjusting  various  parameters  have  on  the  defense 
strategies.  Analysis  of  the  results  obtained  is  presented  in  Chapter  III,  followed 
by  the  conclusion  and  recommendations  for  future  work  that  can  be  done  in  this 
area  of  research  in  Chapter  IV. 
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II.  METHODOLOGY 


A.  SCENARIO 

Consider  a  HVU,  like  an  aircraft  carrier  or  a  supply  container  ship,  moving 
through  waters  with  a  high  number  of  small  ships  in  the  region.  Some  of  these 
small  ships  may  have  no  hostile  intent,  but  hidden  amongst  them  are  a  group  of 
ships  that  are  intending  to  attack  the  HVU.  When  the  HVU  passes  by  their 
location,  these  hostile  ships  can  launch  a  coordinated  attack  against  the  HVU. 
They  maneuver  towards  HVU  at  the  same  time,  and  in  large  numbers.  This  is 
what  we  consider  as  a  swarm  attack.  The  attackers  may  be  equipped  with  small 
caliber  guns  and  some  may  carry  explosive  with  the  intent  of  conducting  a 
suicide  attack.  The  HVU  detects  these  attackers  on  its  radar  and  determines  that 
it  is  under  imminent  attack  based  on  the  number  of  contacts  moving  towards  its 
location  at  high  speed.  The  HVU  then  deploys  multiple  Unmanned  Surface 
Vehicles  (USVs)  equipped  with  weaponry,  to  intercept  and  neutralize  these 
attackers  before  they  can  get  close  enough  to  become  a  threat.  The  number  of 
USVs  (defenders)  deployed  is  small  compared  to  the  number  of  attackers. 

B.  SIMULATION 

The  above  scenario  is  simulated  using  MATLAB®.  The  setup  of  the 
scenarios  to  be  studied  is  to  be  based  upon  real  life  platforms  and  weapons 
systems  to  give  a  realistic  model  and  therefore  a  representative  simulation  of 
actual  hostile  engagement  scenarios.  Current  Fast  Attack  Crafts  (FAC)  and  Fast 
Inshore  Attack  Crafts  (FIAC)  have  max  speeds  in  excess  of  40  knots  (~20  m/s). 
As  such,  the  attackers  are  modeled  to  have  max  speeds  of  45  knots  (~23  m/s), 
while  varying  top  speeds  for  defenders  between  20  knots  (~10m/s)  to  OOknots 
(~30  m/s)  will  be  investigated.  The  hostile  crafts  are  assumed  to  be  armed  with 
small  caliber  guns,  or  carrying  explosives  for  suicidal  attacks.  The  maximum 
effective  ranges  of  these  types  of  weapons  do  not  exceed  1  nmi  (-1.8  km). 
Therefore,  to  be  able  to  intercept  the  hostile  targets  before  they  get  within 
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effective  range  of  their  weapons,  the  defenders  have  to  intercept  these  hostile 
targets  at  least  1nmi  from  the  HVU.  At  the  lowest  defender  speed  we  will  be 
investigating,  the  defenders  will  take  about  3  minutes  to  travel  1  nmi.  During  this 
time,  the  attackers  can  travel  upwards  of  2.25  nmi.  In  other  words,  the  defenders 
have  to  be  deployed  to  intercept  these  hostile  targets  at  least  3.25  nmi  from  the 
HVU.  For  the  simulations,  the  attackers  will  start  at  4.25  nmi  (~8  km)  away  from 
the  HVU  to  provide  a  margin  of  safety.  Modern  sensors  systems  carried  by  Naval 
Fighting  Ships  can  detect  targets  as  far  away  as  100  to  200  km.  So  detecting 
these  targets  would  not  be  a  problem  at  all.  The  weapon  systems  carried  by  the 
defenders  will  be  small  to  medium  caliber  guns  having  max  effective  ranges  from 
Inmi  (-1800  m)  to  about  2  nmi  (-3700  m). 

A  cost  function  is  developed  that  will  allow  us  to  evaluate  the  effectiveness 
of  the  defense  strategies  by  calculating  the  probability  of  the  HVU  surviving  the 
swarm  attack.  We  can  vary  some  of  the  parameters  used  in  the  scenario  to  look 
at  how  some  of  these  parameters  affect  the  survival  rate  of  the  HVU.  Some  of 
the  parameters  we  are  going  to  look  at  are,  guidance  laws  used  by  the  defenders 
to  track  the  attackers,  the  number  of  defenders  and  the  speed  ratio  of  the 
defenders  to  attackers. 

The  following  paragraphs  will  explain  how  the  simulation  is  set  up  to 
simulate  the  motions  of  the  HVU,  attackers  and  defenders,  as  well  as  the 
describing  the  guidance  laws  used  by  the  defenders.  Following  which,  a  matrix  of 
scenarios  is  formed  to  explore  the  effects  of  varying  parameters  in  the  simulation. 

1.  High  Value  Unit  (HVU) 

The  HVU  is  simulated  to  be  a  relatively  slow  moving  vessel  with  limited 
maneuvering  capabilities.  In  the  simulation,  the  HVU  will  start  at  the  origin  (0,  0), 
and  moves  at  a  constant  velocity  of  5  m/s  in  a  straight  line  directly  to  the  north  (in 
the  positive  y  direction).  The  HVU  is  assumed  to  have  detected  the  hostile 
targets  at  a  distance  away  and  is  starting  to  deploy  its  defensive  force  of  USVs. 
The  HVU  will  not  be  taking  evasive  actions  to  avoid  the  attackers  in  this 
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simulation,  since  it  cannot  outmaneuver  the  attackers  anyway.  The  HVU  has  to 
depend  on  the  defensive  force  and  onboard  defensive  capabilities  to  fend  off  the 
attack. 

Figure  1  shows  a  plot  of  the  HVU  motion  in  two-dimensional  space  over 
the  total  simulation  time,  T. 
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2.  Attackers 

The  initial  positions  of  the  40  attackers  are  randomly  generated  with 
uniform  distribution  over  a  rectangular  area  to  the  east  of  the  HVU  (positive  x- 
direction).  This  area  spans  from  x  =  8000  to  x  =  8100,  and  y  =  0  to  y  =  500.  The 
attackers’  trajectories  are  generated  using  a  time-coordinated  path  following 
control  architecture  described  in  Ghabcheloo  et  al,  2009[8],  with  their  speeds 
constrained  to  a  maximum  of  23m/s.  The  trajectories  thus  generated  are 
collision-free,  and  ensures  that  the  attackers  reach  the  HVU  at  the  same  time  for 
a  coordinated  attack.  The  attacker  are  assumed  to  ignore  the  defenders  as  they 
try  to  intercept  them,  as  they  will  be  focused  on  the  task  at  hand  to  destroy  the 
HVU  using  a  suicidal  attack. 
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The  weapons  carried  on  these  attacking  crafts  will  be  limited  to  small  and 
medium  caliber  guns,  as  well  as  explosives,  all  of  which  having  effective  ranges 
of  no  more  than  1800m. 

Figure  2  shows  the  typical  trajectories  of  40  attackers  in  a  coordinate 
attack  against  the  HVU. 


Attackers  trajectories 


Figure  2.  Attackers  trajectories 


3.  Defenders 

The  defenders  are  simulated  to  be  moving  in  formation  alongside  the 
HVU,  and  are  deployed  simultaneously  to  intercept  the  incoming  attackers.  They 
will  start  moving  towards  the  attackers  at  their  maximum  speeds  as  soon  as  they 
are  deployed.  A  divide  and  conquer  tactic  will  be  used,  where  each  defender  will 
be  assigned  a  group  of  attackers  to  intercept  and  neutralize.  The  defenders  will 
use  either  modified  PN  guidance  law  or  pursuit  guidance  law  to  guide  them 
towards  the  group  of  assigned  attackers.  The  guidance  laws  used  will  be 
discussed  in  more  detail  later  in  this  chapter. 

The  field  of  view  (FOV)  for  these  defenders  will  be  ±60°  from  their 

heading.  When  an  attacker  is  within  the  defenders’  FOV  and  weapons’  effective 

range,  the  defenders  will  shoot  at  the  attackers  in  an  attempt  to  neutralize  it.  If 

8 


there  are  multiple  attackers  within  the  defender’s  sight,  the  defender  will  divide  its 
attention  evenly  to  each  of  these  attackers  that  are  in  sight.  This  is  to  ensure  that 
the  defender  can  neutralize  as  many  attackers  as  it  can  before  the  attackers 
maneuvers  past  it. 

Figure  3  shows  a  typical  engagement  scenario  where  five  defenders  try  to 
intercept  40  attackers  that  are  coming  to  attack  the  HVU. 


Defenders  trajectories  to  intercept  attackers 


Figure  3.  Defenders  trajectories  to  intercept  attackers 


4.  Cost  Function 

The  effectiveness  of  the  defense  strategy  is  evaluated  by  estimating  the 
probability  of  the  FIVU  surviving  the  swarm  attack.  A  cost  function,  developed  by 
Claire  Walton  from  UCSC  [9],  derives  the  probability  that  the  HVU  has  survived 
until  time  t  +  At  using  conditional  probabilities.  This  probability  can  be 
represented  as 

p(t  +  At)  =  p(t)  ■  f(t,At)  (2.1) 

where  f{t,At)  is  the  probability  that  the  HVU  survived  in  the  time  interval 
[t,  t  +  At].  By  modeling  the  function  /(t,At),  a  differential  equation  for  p(t)  can 
be  obtained  by  taking  the  limit.  At  ^  0. 
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To  formulate  the  function  /(t,At),  let  us  first  consider  a  single  defender 
protecting  a  HVU  from  a  single  attacker.  Let  p(t)  be  the  probability  that  the  HVU 
survived  until  time  t.  The  goal  is  to  maximize  the  value  of  p(T)  at  the  end  of  the 
simulation  when  t  =  T. 

Let: 

•  q(t)  =  the  probability  that  the  attacker  survived  to  time  t 

•  5a(t)  =  instantaneous  attacker  hit  rate  against  HVU 

•  Sd(t)  =  instantaneous  defender  hit  rate  against  attacker 

The  probability  of  the  attacker  surviving  until  time  t  +  At  can  be  expressed 
as 

q(t  +  At)  =  q(t)(l  -  Sd(t)At)  (2.2) 

This  creates  the  differential  equation 

q(t)  =  -q(t)Sd(t)  (2.3) 


Solving  this  equation  yields 

q(t)  =  (2.4) 

The  effective  attackers  hit  rate  against  the  HVU  at  time  t,  given  that  the 
attacker  survived  until  that  time,  is  therefore  q(t)sa(t).  The  probability  that  the 
HVU  surviving  until  time  t  +  At  can  then  be  expressed  in  a  similar  fashion  to  q(t), 

p(t  +  At)  =  p(t)(l  -  q(t)Sa(t)At)  (2.5) 

which  yields 

p(t)  =  (2.6) 

Substitution  of  q(t)  into  this  expression  and  evaluating  to  time  T,  we  obtain  the 
cost  function, 
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p(T)  =  (2.7) 

To  expand  this  to  a  multi-agent  engagement  scenario,  the  same  method 
can  be  applied,  but  the  hit  rate  of  each  defender  against  each  attacker  has  to  be 
considered.  To  do  this,  we  define  the  following  terms 

•  qi(t)  =  the  probability  that  the  attacker  I  survived  to  time  t 

•  attacker  hit  rate  against  HVU 

•  ^d,/c(t)  =  defender  hit  rate  against  attacker 
Now,  the  probability  that  the  Z-th  attacker  survives  until  time  t  is 

K 

qiit  +  At)  =  qi(t)  j~[(l  -  Sd,fc(t)At)  (2.7) 

k=l 


this  yields 


qj(t)  =  e  ioSLi^d.feWdT 


(2.8) 


Similarly,  we  write  for  the  HVU, 


L 

pit  +  At)  =  p(t)  j~[(l  -  qiit)Sa,iit)At)  (2.9) 

1=1 


This  can  be  made  into  a  differential  equation  if  we  expand  the  product  and 
discard  higher  order  terms  of  At, 


p(t  -I-  At)  =  p(t)  -I- 


Lj 

-^qi(t)' 

1=1 


Saiit)At  +  h.o.t.  —  1 


p(t)  (2.10) 


p(t  -I-  At)  =  p(t)  -I- 


L 

-  ^  qiit)Sa,iit)At  +  h.  0.  t. 

-  1=1 


pit) 


(2.11) 
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The  differential  equation  can  then  be  written 


p(t)  = 


.  k  / 


and  the  solution  evaluated  to  time  T  is 


(2.12) 


T  I 

p(T)  =  g-Jo  (Si=iqi(Osa,i(t))dT  (2.13) 

The  Z-th  attacker  hit  rate  against  the  HVU  is  modeled  using  a  beta  function,  with 
positive  parameters  a,p  >2.  The  function  can  be  connected  smoothly  (C^)  as  a 
piecewise  function, 

=  .0<r(t)<l  (2.14) 

(  0  ,  Q  Is  6 

where  r(t)  is  the  range  of  the  attacker  to  the  HVU.  The  parameters  are  adjusted 
to  reflect  the  hit  rate  of  the  attackers’  weaponry  against  the  HVU.  In  the 
simulations,  the  range  of  this  hit  rate  is  limited  to  800m.  Figure  4  shows  an 
example  of  the  hit  rate  function  modeled  with  the  beta  function. 


50  50 


Figure  4.  Hit  rate  function  modeled  using  beta  function  in  3D 


12 


As  for  the  defenders’  hit  rate  against  the  attackers,  the  weapons  carried  by 
the  defenders  are  assumed  to  be  small  to  medium  caliber  guns.  The  maximum 
effective  ranges  of  such  weapons  are  expected  to  be  around  1800m  to  3700m. 
The  peak  effectiveness  of  the  weapon  is  expected  to  occur  at  a  range  of  around 
100  to  200m  from  the  defender.  This  is  because  if  the  attacker  is  any  closer  to 
the  defender,  the  LOS  rate  to  the  attacker  is  going  to  be  very  high,  especially  for 
a  cross-range  engagement  scenario.  Limitations  of  the  swiveling  rate  of  the  guns 
will  make  it  hard  to  keep  up  with  the  moving  attackers.  The  hit  rate  will  also 
decrease  from  the  peak  value  at  100-200m  at  an  exponential  rate.  This  reflects 
the  fact  that  the  lethality  of  the  shells  fired  at  the  attackers  reduces  with 
increasing  distance  travelled.  Accuracy  of  the  shells  also  decreases  with 
increasing  range  due  to  wind  and  weapon  recoil.  Figure  5  shows  the  expected  hit 
rate  function. 


Defender  Hit  Rate  Function 


Figure  5.  Hit  rate  of  /c-th  defender  against  attacker 

The  hit  rate  function  shown  here  is  modeled  with  a  lognormal  probability 
distribution  curve.  Note  that  the  hit  rate  function  modeled  this  way  does  not  have 
a  finite  range  like  the  beta  functions  have.  The  value  decreases  to  a  very  small, 
but  finite  value  at  very  large  range.  This  leads  to  some  unexpected  results,  which 
will  be  discussed  later  in  the  results  section. 

To  simulate  the  limited  FOV  of  the  weapon  system,  the  hit  rate  function  is 
modified  with  an  angularly  decaying  multiplier,  which  decreases  with  arc  angle. 
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Figure  6  shows  the  multiplier  function  with  such  FOV  limitations  in  three- 
dimensional  space. 
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Figure  6.  Angularly  decaying  rate  function  reflecting  FOV  limitations 

Another  situation  the  hit  rate  function  should  account  for  is  the  fact  that 
when  multiple  attackers  are  within  a  single  defender’s  FOV,  the  defender’s 
attention  is  divided  between  the  attackers.  This  division  of  attention  can  be 
simulated  by  dividing  the  hit  rate  function  by  another  function  that  smoothly 
approximates  the  number  of  attackers  in  its  FOV.  We  can  define 

•  Xfc  =  the  location  of  the  /c-th  defender 

•  xi  =  the  location  of  the  Z-th  attacker 

The  following  sum  of  Gaussian  distributions  can  then  be  used  as  the  division 
function, 

.  /P  -  \\Xk-Xi\\ 

— ? — 

1=1  ^ 

with  the  standard  deviation  a  set  to  a  small  number,  to  approximate  the  number 
of  attackers  within  a  radius  p  of  the  /c-th  defender. 


(2.15) 
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C.  SIMULATION  PARAMETERS 

In  this  section,  we  look  at  the  various  parameters  in  the  simulation  that 
can  influence  the  effectiveness  of  the  defensive  force.  We  will  look  at  three 
parameters,  the  guidance  strategy  the  defenders  use  to  intercept  the  attackers, 
the  number  of  defenders  deployed  and  the  speed  ratio  of  the  defenders’  speed  to 
the  attackers’  speed. 

1.  Guidance  Laws 

Here,  we  are  going  to  look  at  two  guidance  laws.  One  is  based  on  pursuit 
guidance  and  the  other  is  based  on  proportional  navigation  (PN)  guidance. 
Pursuit  guidance  and  PN  guidance  are  common  guidance  laws  used  in  missile 
guidance.  Zarchan  [10]  describes  these  guidance  laws  in  detail.  Modern  missiles 
use  more  complex  forms  of  guidance  laws,  but  in  this  thesis,  we  are  only  going  to 
look  at  the  basic  pursuit  and  PN  guidance  laws  for  a  start.  In  a  nutshell,  pursuit 
guidance  generates  turn  commands  that  points  the  defender  to  the  line-of-sight 
(LOS)  to  the  attacker,  while  PN  guidance  generates  turn  commands  proportional 
to  LOS  rate  to  form  an  intercept  triangle  to  intercept  the  attacker  along  its  path  of 
motion.  Both  of  LOS  angle  and  LOS  rate  can  be  easily  obtained  from  EO  seekers 
that  can  be  equipped  on  a  USV. 

a.  Pursuit  Guidance 

In  the  simulation,  pursuit  guidance  is  implemented  by  having  a 
defender  change  its  heading  to  face  directly  at  the  attacker.  This  is  akin  to  giving 
a  turn  command  that  will  turn  the  defender  an  amount  equal  to  the  LOS  angle  in 
a  single  time  step.  This  simplifies  calculation  and  implementation,  as  we  only 
need  to  find  the  LOS  vector.  When  the  defender  is  up  against  a  swarm  of 
attackers,  the  pursuit  guidance  law  is  not  going  to  work  in  its  basic  form,  as  we 
have  multiple  LOS  vectors  instead  of  just  one.  Therefore,  for  the  defender  to  use 
the  guidance  law,  all  the  LOS  vectors  are  combined  into  one  “effective”  vector 
that  the  defender  will  use.  This  combination  can  be  achieved  by  finding  a 
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weighted  sum  of  the  LOS  vectors.  We  consider  that  in  a  real  life  engagement 
scenario,  priority  should  be  given  to  attackers  that  the  defenders  can  reach  first. 

In  the  simulation,  we  have  the  following  information 


rXwi 

•  Pa  =  =  the  position  of  the  defender 

rXjj  ji 

•  Pa, I  =  [y  J  =  the  position  of  the  Z-th  attacker 


=  the  velocity  of  the  defender 
=  the  velocity  of  the  Z-th  attacker 


•  Va.l  = 


Xa.l 

lya.li 


From  this  we  can  obtain 

•  di,  the  distance  of  the  defender  to  the  Z-th  attacker 

di  =  \\Pd-Pa.i\\  (2.16) 

•  iaj,  the  unit  LOS  vector  of  the  defender  to  the  Z-th  attacker 

ia.i  =  iPd-Pa.i}/di  (2.17) 

•  Vj  is  the  closing  velocity  of  the  defender  to  the  Z-th  attacker 

Vi^ia.i'iVd-Va.i)  (2.18) 


Therefore,  the  unit  velocity  vector  of  the  defender. 


id  = 


^1 


'■a, I 


yL 


I 


=icZ 


a,l 


yL 


(2.19) 


This  places  priority  on  the  attacker  that  the  defender  can  reach  in 
the  shortest  time,  since  —  =  -  where  t,  is  the  time  to  reach  the  Z-th  attacker 

di  ti 

based  on  a  straight-line  intercept  course  of  the  defender  to  the  attacker  along  the 
LOS  vector.  The  lower  the  time,  ti  is,  the  higher  the  weightage  will  be  for  the  Z-th 
attacker.  Figure  7  illustrates  this  method. 
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/ 


Attacker 
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Vl  =  ia/(yd-Va.l) 


yi  ili ;  , 
^1=1 

yi  ^ 

d. 


Figure  7.  Pursuit  guidance  with  weightage  on  LOS  vectors 


Note  that  only  the  attackers  that  are  within  the  FOV  of  the 
defenders,  as  denoted  by  the  shaded  triangle  in  the  figure,  are  included  in  the 
computation. 


b.  PN  Guidance 

For  PN  guidance,  a  turn  command  is  generated  that  is  proportional 
to  the  LOS  rate  of  the  attacker.  This  puts  the  defenders  on  an  intercept  triangle 
where  the  defender  will  intercept  the  attacker  along  its  path  of  motion.  As  we 
have  multiple  attackers  in  this  case,  we  need  to  come  up  with  a  single  turn 
command  for  the  defender  such  that  it  will  track  the  general  swarm  movement, 
while  giving  priority  to  the  closest  threat.  One  possible  way  to  implement  PN 
guidance  for  multiple  attackers  is  to  compute  the  centroid  of  the  attacker  in  the 
defender’s  FOV,  as  well  as  the  effective  velocity  of  that  centroid,  and  apply  PN 
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guidance  to  intercept  this  moving  centroid.  This  centroid  is  to  be  weighted  by  the 
closing  velocity  and  distance  to  each  attacker  to  give  priority  to  attackers  with 
lower  intercept  times.  To  do  this  we  can  find 


The  centroid  of  attackers  in  defender’s  FOV, 


and  the  velocity  of  this  centroid, 


(2.20) 


(2.21) 


From  which,  we  find 


•  da,  the  distance  of  the  defender  to  the  attacker  centroid 

da  =  \\Pd-Pa\\  (2.22) 

•  ia,  the  unit  LOS  vector  of  the  defender  to  the  attacker 
centroid 


'  sin  0, 


cos  0, 


(.Pd  Pa)/da 


Now,  we  can  obtain  the  LOS  rate  by 

(Xd  -  Xa)  sin  Oa  -  (ya  -  fa)  cos  Oa 


0a  = 


The  turn  command  needed,  oj^,  is  then  calculated  by 

0)c  ^0a 

with  the  gain  of  the  PN  guidance  law,  K  =  6 
Figure  8  illustrates  this  method. 


(2.23) 


(2.24) 


(2.25) 
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Figure  8.  Proportional  Navigation  guidance  with  weighted  attacker  centroid 

Again,  note  that  only  attackers  within  the  defender’s  FOV  is 
included  in  the  computation  of  the  required  turn  rate. 

2.  Number  of  Defenders 

By  increasing  the  number  of  defenders  we  have  against  a  swarm,  we  can 
expect  the  defenders  to  take  down  more  attackers  before  they  can  get  to  the 
FIVU.  Flence,  the  chances  of  survival  of  the  FIVU  can  be  expected  to  improve 
with  increasing  number  of  defenders.  Flowever,  we  can  also  expect  that  there 
would  be  some  sort  of  diminishing  return,  where  at  some  point,  increasing  the 
number  of  defenders  is  not  going  to  increase  the  HVU  survival  rate  significantly. 
The  number  of  defenders  deployed  will  be  varied  from  1  to  25. 
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3.  Speed  Ratio 

The  speed  of  the  defenders  will  determine  where  they  can  first  intercept 
the  incoming  attackers.  It  also  affects  the  amount  of  time  that  the  defenders  have 
on  the  target.  We  would  not  want  the  defenders  to  intercept  the  attackers  too 
late.  Otherwise,  they  can  be  close  enough  to  the  HVU  to  pose  a  threat.  However, 
we  would  not  want  the  defenders  to  be  moving  so  fast  that  they  cannot  target  the 
attackers  long  enough  to  take  them  down.  The  speeds  of  the  defenders  will  be 
varied  from  10m/s  to  30  m/s  (about  19knots  to  58knots),  while  the  attackers’ 
speed  will  be  kept  at  23  m/s  (45knots).  This  corresponds  to  speed  ratios  between 
0.4348  and  1.304. 

4.  Test  Matrix 

In  each  of  these  scenarios,  50  randomly  generated  attacker  trajectories 
will  be  evaluated  to  obtain  an  average  performance  of  the  defense  strategy  used. 
This  is  to  make  sure  that  the  performance  of  the  defense  strategy  is  not  a  limited 
to  a  specific  attacker  trajectory.  The  MATLAB®  source  code  written  for  the 
simulation  is  documented  in  Appendix  A  of  this  thesis. 

Table  1  summarizes  the  scenarios  that  are  simulated. 
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Parameter 

Number  of 

defenders 

Number  of 

attackers 

Speed  of 
defenders 
(m/s) 

Speed  of 
attackers 
(m/s) 

Guidance 

Law  used 

1 

40 

25 

23 

Pursuit 

2 

40 

25 

23 

Pursuit 

3 

40 

25 

23 

Pursuit 

4 

40 

25 

23 

Pursuit 

Guidance 

5 

40 

25 

23 

Pursuit 

Law 

1 

40 

25 

23 

PN 

2 

40 

25 

23 

PN 

3 

40 

25 

23 

PN 

4 

40 

25 

23 

PN 

5 

40 

25 

23 

PN 

6 

40 

25 

23 

PN 

7 

40 

25 

23 

PN 

8 

40 

25 

23 

PN 

#  of 

9 

40 

25 

23 

PN 

Defenders 

10 

40 

25 

23 

PN 

15 

40 

25 

23 

PN 

20 

40 

25 

23 

PN 

25 

40 

25 

23 

PN 

5 

40 

10 

23 

PN 

5 

40 

11 

23 

PN 

5 

40 

12 

23 

PN 

5 

40 

13 

23 

PN 

5 

40 

14 

23 

PN 

5 

40 

15 

23 

PN 

5 

40 

16 

23 

PN 

Speed 

5 

40 

17 

23 

PN 

Ratio 

5 

40 

18 

23 

PN 

5 

40 

19 

23 

PN 

5 

40 

20 

23 

PN 

5 

40 

21 

23 

PN 

5 

40 

22 

23 

PN 

5 

40 

23 

23 

PN 

5 

40 

24 

23 

PN 

5 

40 

25 

23 

PN 

5 

40 

30 

23 

PN 
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Parameter 

Number  of 

defenders 

Number  of 

attackers 

Speed  of 
defenders 
(m/s) 

Speed  of 
attackers 
(m/s) 

Guidance 

Law  used 

5 

40 

10 

23 

Pursuit 

5 

40 

11 

23 

Pursuit 

5 

40 

12 

23 

Pursuit 

5 

40 

13 

23 

Pursuit 

5 

40 

14 

23 

Pursuit 

5 

40 

15 

23 

Pursuit 

5 

40 

16 

23 

Pursuit 

Speed 

5 

40 

17 

23 

Pursuit 

Ratio 

5 

40 

18 

23 

Pursuit 

5 

40 

19 

23 

Pursuit 

5 

40 

20 

23 

Pursuit 

5 

40 

21 

23 

Pursuit 

5 

40 

22 

23 

Pursuit 

5 

40 

23 

23 

Pursuit 

5 

40 

24 

23 

Pursuit 

5 

40 

25 

23 

Pursuit 

5 

40 

30 

23 

Pursuit 

Table  1 .  Scenarios  to  be  simulated 
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III.  RESULTS  ANALYSIS 


A.  SIMULATION  RESULTS 

The  measure  of  effectiveness  (MOE)  used  to  evaluate  the  performance  of 
the  defense  strategy  used  is  the  probability  of  the  HVU  surviving  the  swarm 
attack  at  the  end  of  the  simulation,  p(T).  The  larger  the  objective  value  is,  the 
more  effective  the  defense  strategy  is.  The  values  that  p(T)  can  take  are 
between  0  and  1,  with  p(T)  =  0  meaning  that  the  HVU  is  guaranteed  to  be 
destroyed  in  the  swarm  attack  and  p(T)  =  1  meaning  that  the  HVU  is  guaranteed 
to  survive  the  swarm  attack. 

In  this  chapter,  we  look  at  the  performance  of  various  strategies  with 
varying  parameters  as  described  in  Chapter  II  Section  C,  and  analyze  whether 
the  simulation  is  sufficiently  describing  the  engagement  scenario.  In  all  scenarios 
simulated,  there  will  be  40  attackers,  each  having  a  maximum  speed  of  23m/s. 
50  sets  of  initial  attacker  positions  are  generated  randomly  using  uniform 
distribution  over  an  area  from  x  =  8000  to  x  =  8100,  and  y  =  0  to  y  =  500.  This 
corresponds  to  50  different  sets  of  attacker  trajectory.  The  simulation  will  be 
repeated  using  the  same  50  sets  of  attacker  position  for  each  of  the  cases  to  be 
simulated,  and  the  mean  of  the  MOE  obtained  from  the  50  simulations  will  be 
used  as  the  basis  for  comparison  between  cases. 

1.  Guidance  Strategies 

The  performance  of  pursuit  guidance  is  compared  to  that  of  PN  guidance, 
with  the  number  of  defenders  between  1  and  5.  The  speed  of  the  defenders  is 
kept  at  25m/s,  corresponding  to  a  speed  ratio  of  1.1.  The  slight  advantage  in 
speed  for  the  defenders  is  to  allow  the  guidance  laws  to  perform  better  against 
the  attackers.  Table  2  summarizes  the  objective  values  obtained  and  Figure  9 
plots  the  values  for  comparison. 
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Parameter 

Number 

of 

defenders 

Guidance 

Law  used 

MOE 

Mean 

Standard 

deviation 

Guidance  Law 

1 

Pursuit 

0.000314 

0.002160 

2 

Pursuit 

0.092827 

0.041404 

3 

Pursuit 

0.939744 

0.011297 

4 

Pursuit 

0.997371 

0.000407 

5 

Pursuit 

0.999848 

0.000029 

1 

PN 

0.000000 

0.000000 

2 

PN 

0.074510 

0.057751 

3 

PN 

0.171997 

0.282923 

4 

PN 

0.222126 

0.294244 

5 

PN 

0.239250 

0.301488 

Table  2.  Measure  of  Effectiveness  for  pursuit  guidance  and  PN  guidance 
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Figure  9. 


Comparison  of  pursuit  guidance  and  PN  guidance 


The  comparison  between  pursuit  guidance  and  PN  guidance  yields  a 
rather  surprising  result.  Pursuit  guidance,  which  is  usually  considered  inferior  to 
PN  guidance  in  missile  applications,  seems  to  perform  better,  needing  only  three 
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defenders  to  obtain  a  MOE  of  about  0.94.  In  other  words,  the  HVU  survives  94% 
of  the  time  when  three  defenders  are  used  against  a  swarm  of  40  attackers. 
Comparatively,  PN  guidance  can  only  achieve  a  MOE  of  about  0.24  with  five 
defenders,  meaning  the  HVU  survives  only  about  24%  of  the  time.  The  reason 
for  this  result  is  obvious  when  we  look  at  the  trajectory  of  the  defenders.  Figure 
10  shows  a  typical  engagement  scenario  when  pursuit  guidance  is  used. 


Figure  10.  Typical  engagement  scenario  using  pursuit  guidance 

As  can  be  seen  in  the  engagement  scenario,  the  defenders  will  always 
end  up  in  a  tail  chase  against  the  attackers.  This  is  to  be  expected  when  pursuit 
guidance  is  used.  However,  unlike  missile  engagement  scenarios,  the  probability 
of  neutralizing  the  attackers  does  not  depend  on  how  close  the  defenders  can 
get  to  the  attackers,  but  rather  on  the  time  on  target  that  the  defenders  have. 

Once  the  defenders  get  behind  the  attackers,  the  attackers  are  effectively 
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constantly  within  the  defenders’  FOV  and  attacking  range.  Therefore,  the 
probability  of  the  attackers  being  neutralized  increases  greatly. 

Figure  11  shows  the  engagement  scenario  where  five  defenders  are 
deployed  against  a  swarm  of  40  attackers  using  pursuit  guidance.  It  can  be  seen 
from  the  plot  that  the  defenders  were  able  to  neutralize  all  40  attackers  before 
they  can  reach  the  FIVU  by  coming  up  behind  them  and  shooting  them  down. 


Pursuit  Guidance  Engagement  Scenario 


Figure  1 1 .  Pursuit  guidance  engagement  scenario 

PN  guidance,  on  the  other  hand,  gives  rise  to  an  engagement  scenario 
where  the  attackers  are  going  to  pass  in  front  of  the  defenders  and  out  of  their 
FOV  in  a  relatively  short  period.  Since  the  defenders  will  stop  once  no  attackers 
are  within  their  FOV,  they  will  not  give  chase  to  the  attackers  from  behind.  This 
limits  the  amount  of  time  the  defenders  have  in  dealing  with  the  attackers,  and 
subsequently  results  in  a  smaller  chance  of  neutralizing  the  attackers.  A  typical 
engagement  scenario  using  PN  guidance  is  shown  in  Figure  12. 
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Figure  12.  Typical  engagement  scenario  using  PN  guidance 

We  can  see  from  this  engagement  scenario  that  if  attackers  are  able  to 
get  pass  the  defenders,  they  can  reach  the  HVU  without  any  hindrance,  resulting 
in  high  probabilities  of  HVU  being  destroyed.  Figure  13  shows  the  engagement 
scenario  where  five  defenders  are  deployed  against  a  swarm  of  40  attackers 
using  PN  guidance.  The  scenario  shows  that  some  of  the  attackers  were  able  to 
get  past  the  defenders  because  of  the  short  time  period  they  were  passing  into 
the  FOV  of  the  defenders. 
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PN  Guidance  Engagement  Scenario 


Figure  13.  Proportional  Navigation  guidance  engagement  scenario 

2.  Number  of  Defenders 

When  the  number  of  defenders  is  increased,  each  defender  will  have 
fewer  attackers  to  deal  with.  Hence,  the  defenders  should  be  able  to  neutralize 
the  attackers  with  greater  ease.  To  compare  the  effects  of  increasing  the  number 
of  defenders,  we  vary  the  number  of  defenders  from  1  to  25.  The  speed  ratio  is 
kept  at  1.1  for  each  of  the  simulation  run.  Table  3  summarizes  the  objective 
values  obtained  and  Figure  14  plots  the  values  for  comparison. 

The  plot  shows  that  increasing  the  number  of  defenders  does  indeed 
improve  the  chances  that  the  HVU  survives.  One  observation  here  is  that  the 
MOE  increases  almost  proportionally  with  the  increase  of  defenders,  up  to  10 
defenders  when  the  probability  of  HVU  surviving  increases  up  to  around  0.94. 
Increasing  the  number  of  defenders  beyond  10  gives  diminishing  returns,  with 
the  probability  increasing  slightly  to  0.9978  with  15  defenders  (a  99.78%  chance 
of  HVU  surviving).  We  can  deduce  from  this  result  that  each  defender  can  take 
down  four  attackers  efficiently  in  this  particular  engagement  scenario.  Any 
additional  defenders  may  just  mean  that  defenders  are  idling  after  defeating  their 

own  share  of  attackers,  thus  reducing  the  utility  of  these  defenders. 
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Number 

MOE 

Parameter 

of 

defenders 

Mean 

Standard 

deviation 

1 

0.000000 

0.000000 

2 

0.074510 

0.057751 

3 

0.171997 

0.282923 

4 

0.222126 

0.294244 

5 

0.239250 

0.301488 

#  of 

Defenders 

6 

0.434544 

0.334362 

7 

0.497450 

0.314897 

8 

0.698302 

0.213622 

9 

0.827713 

0.134149 

10 

0.935115 

0.065392 

15 

0.997801 

0.001932 

20 

0.999950 

0.000038 

25 

0.999998 

0.000002 

Table  3.  Measure  of  Effectiveness  for  varying  number  of  defenders 
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Figure  14.  Comparison  for  varying  number  of  defenders 
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3.  Speed  Ratio 

The  speed  of  the  defenders  is  now  varied  between  10m/s  and  30m/s.  With 
the  attackers’  speed  at  25m/s,  this  corresponds  to  speed  ratios  between  0.4348 
and  1.3043.  Five  defenders  are  deployed  in  each  case.  Both  pursuit  guidance 
and  PN  guidance  are  evaluated  for  their  effectiveness  at  the  different  speed 
ratios.  Table  4  summarizes  the  objective  values  obtained  and  Figure  15  plots  the 
values  for  comparison. 

The  MOE  when  PN  guidance  is  used  does  not  seem  to  have  a  significant 
difference,  statistically  speaking,  when  the  speed  ratio  is  above  0.74.  However, 
the  defenders  generally  seemed  to  perform  better  at  lower  speed  ratios  between 
0.65  and  0.70.  A  probable  reasoning  behind  this  is  that  at  lower  speed  ratios,  the 
attackers  remain  within  the  defenders’  FOV  for  a  longer  time,  allowing  longer 
dwell  times  for  the  defenders  to  neutralize  the  attackers.  Having  said  that  though, 
the  attackers  are  expected  to  be  intercepted  and  neutralized  at  a  distance  closer 
to  the  HVU  when  the  speed  ratio  is  lower.  This  means  the  safety  margin  is 
smaller  and  the  HVU  may  come  under  fire  if  the  attackers  carry  weapons  with 
longer  range.  At  even  lower  speed  ratios,  we  can  see  that  the  performance 
decreases  dramatically.  One  possible  reason  is  that  at  lower  speeds,  the 
defenders  have  to  aim  at  a  point  far  ahead  of  the  attacker  to  be  able  to  intercept 
the  attackers  using  PN  guidance.  If  they  do  that  though,  the  attackers  may  be  at 
the  edge  of  their  FOV,  lowering  their  hit  rate  against  the  attackers.  The  attackers 
may  be  even  out  of  the  FOV  of  the  attackers  if  the  speed  of  the  defender  is  low 
enough,  and  the  defenders  will  be  unable  to  shoot  at  the  attackers. 

On  the  other  hand,  when  pursuit  guidance  is  used,  a  surprising  result  is 
obtained.  The  MOE  for  simulated  speed  ratios  between  about  0.4348  and  1.304 
when  pursuit  guidance  is  used  is  either  1  or  close  to  one.  This  means  that  the 
HVU  is  almost  guaranteed  to  survive  regardless  of  the  speed  of  the  defenders. 
Additional  simulations  are  run  at  speeds  of  5m/s,  Im/s  and  0.01  m/s, 
corresponding  to  speed  ratio  of  0.2174,  0.0435  and  0.0004  respectively. 
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Guidance 

Law  used 

Speed  of 

Speed 

Ratio 

MOE 

Parameter 

Defenders 

(m/s) 

Mean 

Standard 

Deviation 

10 

0.4348 

0.000 

0.000 

11 

0.4783 

0.000 

0.000 

12 

0.5217 

0.000 

0.000 

13 

0.5652 

0.000 

0.000 

14 

0.6087 

0.020 

0.004 

15 

0.6522 

0.675 

0.326 

16 

0.6957 

0.571 

0.377 

17 

0.7391 

0.230 

0.317 

PN 

18 

0.7826 

0.246 

0.300 

19 

0.8261 

0.168 

0.265 

20 

0.8696 

0.283 

0.295 

21 

0.9130 

0.257 

0.269 

22 

0.9565 

0.344 

0.302 

23 

1.0000 

0.350 

0.316 

24 

1.0435 

0.207 

0.258 

25 

1.0870 

0.239 

0.301 

30 

1.3043 

0.144 

0.230 

Speed 

Ratio 

0.01 

0.0004 

0.188 

0.007 

1 

0.0435 

0.872 

0.004 

5 

0.2174 

1.000 

0.000 

10 

0.4348 

1.000 

0.000 

11 

0.4783 

1.000 

0.000 

12 

0.5217 

1.000 

0.000 

13 

0.5652 

1.000 

0.000 

14 

0.6087 

1.000 

0.000 

15 

0.6522 

1.000 

0.000 

Pursuit 

16 

0.6957 

0.996 

0.000 

17 

0.7391 

0.996 

0.000 

18 

0.7826 

0.999 

0.000 

19 

0.8261 

1.000 

0.000 

20 

0.8696 

1.000 

0.000 

21 

0.9130 

1.000 

0.000 

22 

0.9565 

1.000 

0.000 

23 

1.0000 

1.000 

0.000 

24 

1.0435 

1.000 

0.000 

25 

1.0870 

1.000 

0.000 

30 

1.3043 

0.973 

0.125 

Table  4.  Measure  of  Effectiveness  for  pursuit  guidance  and  PN  guidance  with 

varying  speed  ratios 
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Figure  15.  Comparison  for  varying  speed  ratios  using  pursuit  and  PN  guidance 


This  is  done  to  observe  whether  the  MOE  will  remain  high  when  the 
defenders  are  almost  stationary.  The  MOE  only  started  decreasing  at  a  speed 
ratio  of  0.0435.  Even  at  an  almost  stationary  speed  of  0.01  m/s,  the  MOE  is  at 
0.188.  What  this  means  is  that  even  if  the  defenders  are  barely  moving,  as  long 
as  they  used  pursuit  guidance  against  the  attackers,  the  HVU  can  still  survive 
with  an  18.8%  chance. 

The  reason  that  we  get  such  unrealistic  MOE  for  pursuit  guidance  is  the 
way  the  defenders  hit  rate  against  the  attackers  is  modeled.  As  described  in 
Chapter  II  section  B.4,  the  hit  rate  function  for  the  attackers  is  modeled  with  a 
lognormal  function.  The  way  the  lognormal  function  is  shaped  is  such  that  the 
value  decreases  gradually  from  a  peak  value  in  an  exponential  way.  However, 
the  value  does  not  go  to  zero  even  at  very  large  values.  In  other  words,  there  is  a 
small  but  finite  value  for  the  hit  rate  function  that  is  being  integrated  by  the  part  of 
the  cost  function  that  calculates  the  survival  rate  of  the  attackers.  Since  pursuit 
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guidance  will  always  point  the  defender  in  the  direction  where  the  attackers  are, 
this  small  but  finite  value  is  constantly  causing  the  probability  of  the  attackers 
surviving  to  decrease.  Figure  16  illustrates  this  situation  where  the 

probability  that  the  attackers  survives  is  constantly  decreasing  even  when  they 
should  be  out  of  the  weapon  range  of  the  defenders. 


Plot  of  attackers  changing 
color  as  qj(t)  decreases 


Defendet  Hil  Rale  Function 


Figure  16.  Effect  of  small  but  finite  values  of  hit  rate  function  on  qi(t) 

We  can  see  from  the  trajectory  of  the  attackers  that  the  attackers  almost 
always  remain  within  a  zone  where  the  hit  rate  of  the  defenders  against  them  is 
very  low.  Flowever,  since  the  defenders  are  barely  moving  in  this  case,  the 
closing  speed  of  the  attackers  is  lower  than  when  the  defenders  are  moving  at 
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faster  speeds.  Hence,  the  time  that  the  attackers  are  within  the  defenders’  FOV 
is  considerably  long.  Coupled  with  the  fact  that  the  attackers  are  always  within 
the  defenders’  FOV  due  to  the  nature  of  pursuit  guidance,  the  result  is  that  the 
integration  of  the  low  hit  rate  values  over  a  long  period  of  time  reduces  the 
probability  of  the  attackers  surviving  to  a  point  where  they  are  considered 
neutralized  in  the  simulation. 

This  result  is  not  observed  when  PN  guidance  is  used  because  the 
attackers  will  go  out  of  the  defenders’  FOV  when  the  speed  of  the  defenders  gets 
too  low.  At  higher  speed  ratios,  this  effect  will  not  be  significant  as  the  region  of 
very  low  hit  rate  values  will  be  traversed  in  a  reasonably  short  period  of  time  that 
it  does  not  contribute  much  to  lowering  the  survival  rate  of  the  attackers. 

We  have  to  look  into  altering  the  modeling  of  the  hit  rate  function  of  the 
defenders  to  limit  the  range  at  which  the  function  can  have  a  value.  This  will 
probably  give  more  reasonable  and  more  realistic  results.  This  can  be  done  in 
further  studies  on  this  subject. 
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IV.  CONCLUSION  AND  RECOMMENDED  FUTURE  WORK 


A.  CONCLUSION 

In  this  thesis,  we  have  developed  a  simulation  that  allowed  us  to  simulate 
various  distributed  strategies  that  can  be  used  in  the  defense  of  a  HVU  against  a 
swarm  attack.  A  cost  function  is  formulated  to  calculate  probability  that  the  HVU 
survives  the  swarm  attack.  This  gives  us  a  measure  of  effectiveness  of  each 
strategy  used. 

The  results  have  given  us  great  insights  into  how  some  of  the  defender 
parameters  used  can  influence  the  survival  rate  of  the  HVU.  Pursuit  guidance, 
normally  not  very  effective  in  modern  missile  applications,  has  shown  to  be  a 
rather  effective  strategy,  and  seemingly  even  better  than  the  preferred  PN 
guidance  law. 

Having  larger  number  of  defenders  improves  the  effectiveness  of  the 
defensive  force,  but  beyond  a  certain  number  of  defenders,  the  utility  of  the 
defenders  will  drop.  The  optimal  ratio  of  defenders  to  attackers  is  shown  to  be  1 
to  4  for  the  particular  scenario  that  we  looked  at. 

The  simulations  also  showed  that  having  high  defender  speeds  might 
hinder,  not  help,  the  defense  operation.  A  moderate  speed  that  allows  sufficient 
time  on  target,  but  yet  able  to  intercept  attackers  at  a  reasonable  distance  from 
the  HVU  should  be  chosen. 

The  results  where  pursuit  guidance  at  low  speeds  perform  surprisingly 
well  reveals  that  the  hit  rate  function  used  in  the  simulations  still  needs  some 
more  considerations  to  make  the  analysis  more  accurate. 

B.  RECOMMENDED  FUTURE  WORK 

1.  Other  Guidance  Methods  or  Strategies 

The  simulation  developed  in  this  thesis  provides  a  framework  in  which 
different  guidance  methods  or  defensive  strategies  can  be  evaluated  for  their 
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effectiveness.  For  example,  Shin,  2011  [6]  describes  a  Earliest  Intercept 
Geometry  Guidance  Law  (EIGGL)  that  can  be  used  for  target  allocation  to 
minimize  the  distance  that  the  attacker  can  cover,  and  hence  maximizes  the 
distance  of  the  attackers  to  the  HVU. 

An  impact  angle  control  guidance  law  is  also  described  in  the  same  paper. 
This  guidance  law  optimizes  the  impact  angle  between  the  defender  and  the 
attacker  during  the  time  of  intercept,  effectively  setting  them  up  in  a  head-on 
collision  course,  which  increases  the  chance  of  interception.  These  laws  can  be 
evaluated  for  their  effectiveness  in  the  swarm  attack  scenario  that  we  have  set 
up. 


2.  Different  Engagement  Scenarios 

In  this  thesis,  we  have  only  looked  at  a  particular  engagement  scenario 
where  a  HVU  is  moving  in  a  straight  line  at  a  relatively  slow  speed,  with  attackers 
moving  aggressively  towards  it  from  one  direction.  A  HVU  can  come  under  a 
swarm  attack  in  many  other  possible  situations.  For  example,  attackers  can 
come  from  multiple  directions,  and  may  maneuver  to  avoid  the  defenders  that  are 
deployed  against  them.  The  HVU  can  be  stationary,  like  an  oilrig  or  even  an  on¬ 
shore  facility.  All  these  scenarios  can  be  modeled  and  simulated  within  this 
framework,  and  different  defensive  strategies  can  be  evaluated  to  find  out  their 
effectiveness  in  such  scenarios. 

3.  Intent  Recognition 

At  the  initial  stage  of  the  research  for  this  thesis,  a  real  time  graphical 

simulation  with  intent  recognition  algorithm  was  being  developed  in  the  University 

of  Nevada,  Reno  (UNR).  The  intention  was  to  integrate  this  simulation  with  the 

simulation  developed  here  to  simulate  the  engagement  in  real  time  and  display 

the  engagement  in  3D  graphics.  More  importantly,  the  intent  recognition 

algorithm  would  allow  us  to  study  the  reaction  of  the  defensive  force  when  the 

hostile  targets  are  hidden  among  neutral  ones.  Daniel  Bigelow  from  UNR,  came 

down  to  NPS  to  install  the  simulation  for  the  integration.  Unfortunately,  even 
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though  the  simulation  was  successfully  installed,  the  simulation  was  unstable  due 
to  the  different  hardware  used  here  in  NPS.  Due  to  time  constraints,  this  part  of 
the  thesis  had  to  be  abandoned,  but  future  work  can  be  done  to  further  explore 
the  possibilities  of  integrating  the  simulations. 

4.  Hit  Rate  Function  Modeling 

As  discussed  in  the  results  section,  the  hit  rate  functions  used  currently 
may  give  us  MOE  that  are  misleading  due  to  the  small  finite  values  of  the  hit  rate 
function  at  large  ranges.  Further  work  is  needed  to  come  up  with  a  hit  rate 
function  that  has  a  finite  range  value  and  describes  range  effectiveness  of  real 
life  weapons  more  accurately. 
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APPENDIX  A.  MATLAB®  SIMULATION  CODE 


This  appendix  contains  the  list  of  script  files  and  function  codes  that  were 
used  in  the  simulation. 

A.  SIMULATION  RUN  SCRIPT  FILES 


g, _ 

o 

%  File:  Evaluate_Scenario_Performance .m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Sets  up  the  parameters  for  cost  evaluation  algorithm 
%  -  Runs  trajectory  generation  scripts 

%  -  Computes  cost  and  plots  attacker  and  defender  trajectories 
%  Inputs:  Discretization  time  and  space,  method  of  evaluation,  number 
%  of  attackers,  number  of  defenders 

%  Outputs:  Objective  value 

g, _ _ 

o 

function  ObjValue  =  Evaluate  Scenario  Performance (Discretization, 

Methods , 

Numberof Attackers , 

Numberof Defenders ) 


global  CONSTANTS  OFFLINE_TRAJECTORIES 
PDF_VALUES  MESHED_PDF_VALUES 

DISCRETIZATION_VALUES  MESHED_DISCRETIZATION_VALUES 
DIFFERENTIATION_MATRICES 

INTEGRATION_WEIGHTS  MESHED_INTEGRATION_WEIGHTS ; 
global  pO  N  %  To  interface  with  Kaminer  code  without  editing 
warning ( ' off ' ,  ' all ' ) 

%  This  suppresses  the  warning  that  future  versions  of  MATLAB®  will  not 
%  support  evaluating  scripts  with  feval,  which  is  what  I'm  using  to 
%  evaluate  arbitrary  versions  of  Ding's  heuristic 

addpath ( ' . /performance_related_f unctions ' ,  ' -begin ' ) ; 

addpath ( ' . /performance_related_functions/Parameterized  Control_Kernel ' ) 
addpath ( ' . /perf ormance_related  functions/hit_rate_functions ' ) 
addpath ( ' ./performance  related_functions/cline ' ) 
addpath ('. /code ' ,  '-begin')  %Ding's  code 

Heuristic  =  ' PlotOneManyBadGuysVl 1 ' ; 

%  PDF  choices 

CONSTANTS . ParameterSpace . PFD  Choices 

= { ' Independent ' , ' Uniform ' , ' Uniform ' } ; 


% - Simulation  Time  Range - % 

CONSTANTS. Time. T0=0;  %start  time 
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% - Parameter  Ranges - % 

g,  _  ___  __  __  __  __  _ 

o 

%  wo  and  WF  should  be  IxParameterSpace . Dimension 
%  arrays  with  the  starting  points  and  ending  points 
%  respectively  of  each  parameter's  domain. 

%  When  attackers  have  identical  parameter  ranges, 

%  the  shortcut  of  listing  just  the  number  of 
%  parameters  per  attacker  can  be  taken. 

g, _ 

o 


g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 


CONSTANTS . ParameterSpace . W0= [ 80 00 , 0 ] ;  %start  of  each  parameter 

CONSTANTS . ParameterSpace . WF= [ 81 00 , 500 ] ;  %end  of  each  parameter 


9-__  ___  ___  _  _  __  _  _____  __  _ 

o 

%  The  following  variables  are  not  meant  to  be  modified. 

%  They  are  set  in  the  function  call  and  just  relabeled  here  to  match 
%  with  previous  code. 

g,  __  _  _  _  _  _  ___  _  _  ___  _  _  _ 

o 

N  =  Discretization ( 1 ) ; 

CONSTANTS. N  =  N;  %  This  redundant  declaration  is  here,  because  the 
%  heuristics  and  Kaminer  code  need  N  as  a  global  variable  to  be  run 
%  with  minimal  editing  (to  be  run  without  rewriting  them  as  function 
%  calls),  but  the  control  code  and  many  subroutines  use  CONSTANTS. N 

%  Number  of  probabilistic  parameters  per  attacker.  In  this  case, 

%  starting  position  (in  two  dimensions) . 

CONSTANTS . ParameterSpace . Dimension  =2; 

%  Number  of  attackers 

%  This  is  part  of  the  substruct  of  attacker  specific  constants,  for 
%  historical  reasons. 

CONSTANTS .ATTACKERS .Na  =  Numberof Attackers ; 

%  Number  of  searchers 
CONSTANTS. Ns  =  Numberof Defenders ; 

g,  ____  ____  __  __  _  _ 

o 


g, 

o 


%  Set  the  constants  for  the  hit  rate  functions  and  the  filenames 


Calibrate  Hit  Rate  Functions 


—  —  —  —  —  —  —  —  —  —  —  —  —  —  —  _____  ___  _ 

o 

%  Create  DIFFERENTIATION_MATRICES ,  INTEGRATION_WEIGHTS , 

%  MESHED_INTEGRATION_WEIGHTS,  DISCRETIZATION  VALUES,  and 
%  MESHED  DISCRETIZATION  VALUES. 


g, 

o 


U.-Lo^  V  oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

disp (Heuristic) ; 

disp(' - '); 

disp ( 'Discretization :  '  )  ; 
disp (Discretization)  ; 

Calculate  Methods (Discretization, 
disp(' - 


Methods ) 


'  ) 


g, 

o 
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[x_a_i,  x_d_a_i,  x_0_a_i]  = 

RunDingsSimulation (Heuristic, Numberof Defenders , 
Numberof Attackers ) ; 

0, _ 

o 

for  s  =  1 : Numberof Attackers 

X  =  interpl (x_a_i { s }  ( : ,  3 )  ,  [x_a_i { s }  ( : , 1 ) , x_a_i { s }  ( : , 2 ) ] , 

DISCRETIZATI0N_VALUES{1}, 'spline' ) ; 
x_a_i{s} ( : , 3)  =  DISCRETIZATION_VALUES { 1 } ( : ) ; 
x_a_i { s } ( : , 1 )  =  X ( : , 1 ) ; 
x_a_i { s }  ( : , 2 )  =  X ( : , 2 ) ; 

end 

for  s  =  1 : Numberof Defenders 

X  =  interpl (x_d_a_i { s }  ( : ,  3 )  ,  [x_d_a_i { s }  ( : , 1 ) , x_d_a_i { s }  ( : , 2 ) ] , 

DISCRETIZATI0N_VALUES{1}, 'spline' ) ; 
x_d_a_i{s} ( : , 3)  =  DISCRETIZATION_VALUES { 1 } ( : ) ; 
x_d_a_i { s } ( : , 1 )  =  X ( : , 1 ) ; 
x_d_a_i { s }  ( : ,  2 )  =  X ( : , 2 ) ; 

end 

X  =  interpl (x_0_a_i ( : , 3 ) , [x_0_a_i ( : , 1 ) , x_0_a_i ( : , 2 ) ] , 
DISCRETIZATI0N_VALUES{1}, 'spline' ) ; 
x_0_a_i(:,3)  =  DISCRETIZATION_VALUES { 1 } ( : ) ; 
x_0_a_i ( : , 1 )  =  X ( : , 1 ) ; 
x_0_a_i ( : , 2 )  =  X ( : , 2 ) ; 

9~—  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  — 

o 


%  Plot  various  desired  plots. 

q, _ _ _ 

o 


Result_Plots ; 

q,  _ _ 

o 


%  Calculate  the  probability  function  just  for  the  a  i  trajectories 


[p,  q]  =  Probability  Function (x  a  i , x  d  a  i ,  x  0  a  i )  ; 

%  Probability  that  the  HVU  does  not  survive  at  time  T. 
ObjValue  =  l-p(end); 

q, _ _ _ _ _ _ _ _ _ 

o 
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_________  _ 

o 

%  File:  Run  Scenarios. m 

%  Compiler :~MATLAB®  v7 . 10 . 0 . 499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Runs  script  to  load  50  pre-generated  attackers'  position 
%  -  Runs  Evaluate_Scenario  Performance . m  with  appropriate  parameters. 

%  -  Logs  evaluated  performance  for  the  50  runs 
%  Inputs:  None 

%  Outputs:  Array  of  objective  values 

g, _ _ _ _ _ _ _ _ _ _ _ 

o 

global  pO 

ObjValues  =  zeros ( 50 , 1 ) ; 
for  sen  =  1:50 

pO  =  load ([' Scenario '  int2str(scn)  ' . dat ' ] ) ; 

ObjValues (sen)  =  Evaluate_Scenario_Perf ormance ( [500, 1,1], [0,0,0], 
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B.  ATTACKER  &  DEFENDER  TRAJECTORY  GENERATION  CODE 


g, _ _ 

o 

%  File:  RunMoreThanOneBadGuy . m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Generates  optimized  attacker  trajectory  based  on  given 
%  limitations 

%  Inputs:  None 

%  Outputs:  Time  for  optimized  trajectory 

9-  __  _  _______ 

o 

global  vShip  psiShip  vO  vf  vmin  vmax  psidotmax  pO  al  lambdaO  lambdaf  tf 
pOship  BadGuyTotal  N  DefenderTotal 
global  p  swarm  Defender  p  ship  Dindex  tf 

%  Initial  conditions  for  the  ship 
vShip  =  5;  %m/sec 

psiShip  =  pi/2;  %rad 

pOship  =  [ 0 ; 0 ]  ; 

%  Initial  conditions  for  the  bad  guys 

%  Initial  and  final  velocities  of  the  bad  guys  in  m/sec 
vO  =  23;  vf  =  23; 

%  Limits  for  the  bad  guys 
%  Max  and  min  speeds  of  the  bad  guy 
vmin  =  1;  vmax  =  23; 

%  Turn  rate  limit  for  the  bad  guy 
psidotmax  =  0.1; 

%  Initial  guess  on  the  hit  time 
tf  =  norm(p0(:,l)  -  pOship) /vmax; 
for  i  =  2 : BadGuyTotal 

tf  =  min (tf, norm (pO (:, i)  -  pOship) /vmax) ; 

end 

options  =  optimset ( ' TolFun ' , . 1 , ' maxiter ' , 1 00 , ' MaxFunEvals ' , 1 00 , 

' Display ' ,  'on'); 

[x, fval, exitflag, output]  = 

fminsearch (@MoreThanOneBadGuyCost, tf , options) ; 

%  This  ends  the  initialization  of  BadGuy  part 
tf  =  X (1)  ; 
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g,_ _ _ _  _ _ _ _  _  _ _ _  _ _ _  _ _ _ _ _ _ 

o 

%  File:  MoreThanOneBadGuy .m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Computes  cost  for  a  single  attacker 
%  Inputs:  Number  of  attackers,  final  time 

%  Outputs:  Cost  for  single  attacker  trajectory 

0,  __  _____  ____  ____ 

o 

function  J  =  MoreThanOneBadGuy (BadCuyNumber, tfO ) 

global  al  vShip  psiShip  vO  vf  vmin  vmax  psidotmax  pO  pOship  N 

%  Time  step 
dt  =  tfO/N; 

%  Extract  optimization  parameters 
xO  =  pO ( : , BadCuyNumber) ; 

psiO  =  atan2 (xO (2) -pOship (2) , xO (1) -pOship (1) ) -pi; 

%  Initialize 

xf  =  pOship  +  [vShip*cos (psiShip) ; 

vShip*sin (psiShip) ]* (tfO- (BadCuyNumber  -  l)*dt)'; 
xpO  =  [cos(psiO);  sin(psiO)]; 

perror  =  [cos (psiShip)  sin (psiShip) ; -sin (psiShip)  cos (psiShip) ] 

* (xO  -  pOship) ; 

if  perror (2)  >=  0 

psif  =  psiShip  -  pi/2; 

else 

psif  =  psiShip  +  pi/2; 

end 

xpf  =  [cos (psif ) ;  sin  (psif)]; 

xppO  =  zeros (2,1);  xppf  =  zeros (2,1); 

xpppO  =  zeros (2,1);  xpppf  =  zeros (2,1); 

%  Define  speed  profile 
lambdaO  =  vO/norm (xpO ) ; 
lambdaf  =  vf /norm (xpf ) ; 

%  Evaluate  UAV  path  coefficients 
if  (abs (lambdaO  -  lambdaf)  <  le-6) 
tauf  =  Iambda0*tf0; 

else 

tauf  =  (lambdaf  -  lambdaO) *tfO/log (lambdaf/lambdaO) ; 

end 

%  The  7th  order  coefficients  below  are  obtained  using  the  same  steps  as 
%  in  Ghabcheloo  paper 
%A  =  [10000000; 

%  01000000; 

%  00200000; 

%  00060000; 

%  1  tauf  tauf''2  tauf^3  tauf''!  tauf^S  tauf''6  tauf^7; 

%  0  1  2*tauf  3*tauf^2  4*tauf^3  5*tauf^4  6*tauf''5  7*tauf^6; 


44 


g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 


002  6*tauf  12*tauf^2  20*tauf^3  30*tauf"4  42*tauf^5; 
0  0  0  6  24*tauf  60*tauf^2  120*tauf-'3  210*tauf^4]  ; 

b  =  [xO  xpO  xppO  xpppO  xf  xpf  xppf  xpppf ] ' ; 

al  =  inv (A) *b; 


al  =  [xO  . . . 

xpO  . .  . 

xppO/2  .  .  . 
xpppO/6  . . . 

(35*xf) /tauf^4  -  (35*x0) /tauf-"4  -  (20*xp0) /tauf'"3  - 

( 15*xpf ) /tauf ^3  -  (5*xpp0) /tauf^2  +  (5*xppf ) / (2*tauf^2) 
(2*xppp0) / (3*tauf)  -  xpppf /( 6*tauf)  ... 

(84*x0) /tauf^S  -  (84*xf) /tauf^S  +  (45*xp0) /tauf"4  + 

(39*xpf) /tauf^4  +  (10*xpp0) /tauf^3  -  (7*xppf) /tauf^3  + 
xpppO/tauf ^'2  +  xpppf/  (2*tauf''2)  .  .  . 

(70*xf) /tauf^e  -  (70*x0) /tauf"6  -  (36*xp0) /tauf^S  - 
(34*xpf) /tauf'^S  -  (15*xpp0)  /  (2*tauf''4)  + 

( 13*xppf ) / (2*tauf ^4 )  -  (2*xppp0) / (3*tauf^3)  - 
xpppf/ (2*tauf"3)  ... 

(20*x0) /tauf^7  -  (20*xf) /tauf''7  +  (10*xp0) /tauf'^S  + 

(10*xpf) /tauf^6  +  (2*xpp0) /tauf ^5  -  (2*xppf ) /tauf ^5  + 
xpppO/ ( 6*tauf "'4 )  +  xpppf /( 6*tauf"'4 )]  ; 

V  max  =  max(vO,vf);  v  min  =  min(vO,vf);  psidot  max  =  0; 

for  n  =  0:N 

t  =  n*dt; 

if  (abs(lambdaO  -  lambdaf)  <  le-6) 
tau  =  t/tfO*tauf; 

else 

tau  =  tauf*  (  (lambdaf/lambdaO) (t/tfO)  -  1)  *lambdaO/ (lambdaf- 

lambdaO ) 

end 

%  Compute  UAV  path  and  its  derivatives 
p  =  zeros (2,1); 
for  i  =  1:8 

p  =  p  +  al  ( : ,  i)  *tau'' (i-1)  ; 

end 

p_p  =  zeros (2,1); 
for  i  =  2 : 8 

p_p  =  p_p  +  (i-1) *al ( : , i) *tau^ (i-2) ; 

end 

p_pp  =  zeros (2,1); 
for  i  =  3 : 8 

p_pp  =  p_pp  +  ( i-1 )  *  ( i-2 )  *al  ( : ,  i )  *tau'' ( i-3 )  ; 

end 

%  Compute  speed 

lambda  =  lambdaO  +  ( lambdaf-lambdaO ) *tau/tauf ; 

V  =  lambda*norm (p  p) ; 
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%  Compute  signed  curvature  in  2D 

curv  =  (p_pp (2) *p_p (1)  -  p_pp (1) *p_p (2) ) /norm(p_p) ^3; 
psidot  max  =  max(psidot  max, v*curv) ; 

V  max  =  max(v,v  max); 

V  min  =  min(v,v  max); 

end 

J  =  tfO  +  10* (psidot  max  -  psidotmax) ^2 /psidotmax^2  +  (v  min  - 
vmin) '^2 /vmin"'2  +  le6*  (v  max  -  vmax) '^2 /vmax''2  ; 

return 


g,  _ _ _  _ _ _  _ _ _  _ _ _  _ _ _  _ _ _  _ 

o 

%  File:  MoreThanOneBadGuyCost  .m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Computes  total  cost  for  all  attackers 
%  Inputs:  Final  time 

%  Outputs:  Total  cost  for  all  attackers 

o  _  _  ____  ___ 

o 

function  J  =  MoreThanOneBadGuyCost (tfO) 
global  BadGuyTotal 
J  =  t  f  0 ; 

for  i  =  1 : 1 : BadGuyTotal 

J  =  J  +  MoreThanOneBadGuy (i,  tf 0)  ; 

end 

return 
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g,_ _ _ _ _ _ _ _ _  _  _ _ _ _ _ _ _ _ 

o 

%  File:  MoreThanOneBadGuyData  .m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Calculate  attacker  position  and  velocity 
%  Inputs:  Number  of  attackers,  final  time,  current  simulation  time 
%  Outputs:  Position,  speed,  turn  rate,  velocity  of  attacker 

9-  _________________ 

o 

function  [p, v, psidot, vel]  =  MoreThanOneBadGuyData (BadCuyNumber, tfO, t) 

global  vShip  psiShip  vO  vf  pO  pOship  N  %  vmin  vmax  psidotmax 

%  Time  step 
dt  =  tfO/N; 

%  Extract  optimization  parameters 
xO  =  pO ( : , BadCuyNumber) ; 

psiO  =  atan2 (xO (2) -pOship (2) , xO (1) -pOship (1) ) -pi; 

%  Initialize 

xf  =  pOship  +  [vShip*cos (psiShip) ;  vShip*sin (psiShip) ] * (tf 0- 

(BadCuyNumber  -  l)*dt) ' ; 

xpO  =  [cos(psiO);  sin(psiO)]; 

perror  =  [cos (psiShip)  sin (psiShip) ; -sin (psiShip)  cos (psiShip) ]* (xO- 

pOship) ; 


if  perror (2)  >=  0 

psif  =  psiShip  -  pi/2; 

else 

psif  =  psiShip  +  pi/2; 

end 

xpf  =  [cos (psif ) ;  sin  (psif)]; 

xppO  =  zeros (2,1);  xppf  =  zeros (2,1); 

xpppO  =  zeros (2,1);  xpppf  =  zeros (2,1); 

%  Define  speed  profile 
lambdaO  =  vO/norm (xpO ) ; 
lambdaf  =  vf /norm (xpf ) ; 

%  Evaluate  UAV  path  coefficients 
if  (abs (lambdaO  -  lambdaf)  <  le-6) 
tauf  =  Iambda0*tf0; 

else 

tauf  =  (lambdaf  -  lambdaO) *tf0/log (lambdaf/lambdaO) ; 

end 

%  The  7th  order  coefficients  below  are  obtained  using  the  same  steps  as 
%  in  Ghabcheloo  paper 
%A  =  [10000000; 

%  01000000; 

%  00200000; 

%  00060000; 

%  1  tauf  tauf''2  tauf^3  tauf''!  tauf^S  tauf''6  tauf^7; 

%  0  1  2*tauf  3*tauf^2  4*tauf^3  5*tauf^4  6*tauf''5  7*tauf^6; 
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g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 


002  6*tauf  12*tauf^2  20*tauf^3  30*tauf"4  42*tauf^5; 
0  0  0  6  24*tauf  60*tauf^2  120*tauf-'3  210*tauf^4]  ; 

b  =  [xO  xpO  xppO  xpppO  xf  xpf  xppf  xpppf ] ' ; 

al  =  inv (A) *b; 


al  =  [xO  . . . 

xpO  . .  . 

xppO/2  .  .  . 
xpppO/6  . . . 

(35*xf) /tauf^4  -  (35*x0) /tauf-"4  -  (20*xp0) /tauf'"3  - 

( 15*xpf ) /tauf ^3  -  (5*xpp0) /tauf^2  +  (5*xppf ) / (2*tauf^2)  - 
(2*xppp0) / (3*tauf)  -  xpppf /( 6*tauf)  ... 

(84*x0) /tauf^S  -  (84*xf) /tauf^S  +  (45*xp0) /tauf"4  + 

(39*xpf) /tauf^4  +  (10*xpp0) /tauf^3  -  (7*xppf) /tauf^3  + 
xpppO/tauf ^'2  +  xpppf/  (2*tauf''2)  .  .  . 

(70*xf) /tauf^e  -  (70*x0) /tauf"6  -  (36*xp0) /tauf^S  - 
(34*xpf) /tauf'^S  -  (15*xpp0)  /  (2*tauf''4)  + 

( 13*xppf ) / (2*tauf ^4 )  -  (2*xppp0) / (3*tauf^3)  - 
xpppf/ (2*tauf"3)  ... 

(20*x0) /tauf^7  -  (20*xf) /tauf''7  +  (10*xp0) /tauf'^S  + 

(10*xpf) /tauf^6  +  (2*xpp0) /tauf ^5  -  (2*xppf ) /tauf ^5  + 
xpppO/ ( 6*tauf "'4 )  +  xpppf /( 6*tauf"'4 )]  ; 

if  (abs (lambdaO  -  lambdaf)  <  le-6) 
tau  =  t/tfO*tauf; 

else 

tau  =  tauf*  (  (lambdaf/lambdaO) (t/tfO)  -  1)  *lambdaO/ (lambdaf- 
lambdaO) ; 
end 

%  Compute  UAV  path  and  its  derivatives 
p  =  zeros (2,1); 
for  i  =  1 : 8 

p  =  p  +  al  ( : ,  i)  *tau'' (i-1)  ; 

end 

p_p  =  zeros (2,1); 
for  i  =  2 : 8 

p_p  =  p_p  +  (i-1) *al ( : , i) *tau^ (i-2) ; 

end 

p_pp  =  zeros (2,1); 
for  i  =  3 : 8 

p_pp  =  p_pp  +  (i-1)  *  (i-2)  *al  ( : ,  i)  *tau'' (i-3)  ; 

end 

%  Compute  speed 

lambda  =  lambdaO  +  ( lambdaf-lambdaO ) *tau/tauf ; 

V  =  lambda*norm (p  p) ; 
vel  =  lambda*p  p; 

%  Compute  signed  curvature  in  2D 

curv  =  (p_pp (2) *p_p (1)  -  p_pp (1) *p_p (2) ) /norm(p_p) ^3;  %  THIS  IS  IN  M 

psidot  =  v*curv; 

end 
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g,_ _ _ _ _  _ _ _  _ _ _ _ _  _ _ _  _ _ _ _ _ _ 

o 

%  File:  PlotOneManyBadGuysVll .m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  Each  defender  calculates  a  weighted  centroid  and  effective  velocity 
%  -  for  the  attackers  in  its  FOV. 

%  -  A  turn  rate  command  is  generated  using  PN  guidance  using  the 
%  -  calculated  centroid  and  velocity. 

%  -  This  version  divides  up  the  attackers  and  assigns  them  to  each 
%  -  defender.  The  defenders  will  stop  once  they  lose  sight  of  attackers 
%  -  assigned  to  them. 

%  -  This  version  also  includes  a  neutralizing  algorithm,  in  which  an 
%  -  attacker  is  neutralized  when  it  stays  within  the  FOV  and  attacking 
%  -  range  of  a  defender  for  a  set  amount  of  time. 

%  Inputs:  Initial  ship  and  swarm  parameters 
%  Outputs:  Defender  trajectory 

o  _______  _  __  __ 

o 

global  vShip  psiShip  vO  vf  N  pO  DefenderTotal 
global  p  swarm  Defender  p  ship  Dindex  tf 

%  PlotOneManyBadGuysVll  -  Guidance  Version  11 

dt  =  tf/N;  %  Time  Step 
p_ship  =  zeros (N+1 , 3 ) ; 
for  i  =  0 : N 
t  =  i*dt; 

p_ship ( i+1 , : )  =  [t  pOship'  +  [vShip*cos (psiShip) 
vShip*sin (psiShip) ] *t] ; 
end 

%  Assigning  attackers  to  defenders 
%  This  algorithm  assumes  attackers  >=  defenders 

%  Attacker  groups  are  assigned  to  defenders  based  on  their  position 
Assign=zeros (DefenderTotal,  2)  ; 

Assign (DefenderTotal,  1) =1; 

Assign (DefenderTotal, 2) =floor (BadGuyTotal/DefenderTotal) ; 
for  i=DefenderTotal-l : -1 : 1 

Assign ( i , 1 ) =Assign ( i+1 , 2 ) +1 ; 

Assign (i, 2) =Assign (i  +  1, 2) +floor (BadGuyTotal/DefenderTotal)  ; 

end 

for  i=l :mod (BadGuyTotal , DefenderTotal) 

Assign ( 1 : i , 2 ) =Assign ( 1 : i , 2 ) +1 ; 

end 

for  i=l :mod (BadGuyTotal, DefenderTotal) -1 
Assign ( 1 : i , 1 ) =Assign ( 1 : i , 1 ) +1 ; 

end 

%  Centroid  of  initial  position  of  attackers 
APos  =  zeros (2, DefenderTotal) ; 
for  j =1 : DefenderTotal 

for  i  =  Assign (j , 1) :Assign (j , 2) 

APos(:,j)  =  APos(:,j)  +  p0(:,i); 

end 

APos ( : , j )  =  APos (:, j )/ (Assign (j , 2) -Assign (j , 1) +1) ; 

end 
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%  Defenders  parameters 

MaxOmega  =1;  %  Max  turn  rate 

Dfov  =  120;  %  Field  of  view  (deg) 

Dv  =  25;  %  Speed  (m/s) 

%  Initialize  matrices 

range  max  =  0;  v  max  =  max(v0,vf);  v  min  =  min(v0,vf);  psidot  max  = 
p_test  =  zeros (N+1 , 2 ) ; 

p  swarm  =  zeros (N+1, (BadGuyTotal*2 ) +1 ) ;  v  swarm  = 
zeros (N+1 ,  BadGuyTotal  +  1 )  ;  psi_dot  =  zeros (N+1 , BadGuyTotal  +  1 ) ; 
pcombined  =  zeros ( 1 ,  BadGuyTotal*2 ) ;  vcombined  =  zeros ( 1 , BadGuyTotal ) 
psidotcombined  =  zeros ( 1 , BadGuyTotal ) ; 

Dpos  =  zeros (2, DefenderTotal) ;  Dpsi  =  zeros (DefenderTotal, 1) ; 

Di  =  zeros (2 , DefenderTotal ); Defender  =  zeros (N+1,  DefenderTotal*2  +  l ) 
VelD  =  zeros (2, DefenderTotal) ; 

Range  =  zeros (DefenderTotal , N+1 , BadGuyTotal+1 ) ;  Dindex  = 
zeros (DefenderTotal,  1)  ; 

Alog  =  zeros (N+1,  DefenderTotal*2  +  l )  ;  HeadingLog  =  zeros (N+1, 
DefenderTotal+1 ) ; 

NeutAtt  =  zeros (BadGuyTotal ,  DefenderTotal  +  1); 

SimEnd  =  0; 

%  Initial  positions  and  heading  for  defenders 
for  i  =  1 : DefenderTotal 

Dpos(:,i)  =  pOship  +  (i-l)*[0;25]  +  [10;0]; 

%  Velocity  vector  pointed  towards  centroid  of  attackers 

VelD ( : , i )  =  Dv* (APos ( : , i ) -Dpos ( : , i ) ) /norm (APos ( : , i ) -Dpos ( :  ,  i ) ) ; 

Dpsi(i)  =  pi/2  -  atan2 (VelD (2 , i) ,  VelD(l,i)); 

end 

for  i  =  0 : N 
t  =  i*dt; 

p_ship ( i+1 , : )  =  [t  pOship'  +  [vShip*cos (psiShip) 

vShip*sin (psiShip) ] *t] ; 


WPos  =  zeros (2,  DefenderTotal); 

WVel  =  zeros (2,  DefenderTotal); 

W  =  zeros (DefenderTotal, 1) ; 

for  j  =  1 : BadGuyTotal 

[p, V, psidot, vel]  =  MoreThanOneBadGuyData ( j  ,  tf ,  t)  ; 
pcombined (1,  ((2*j)-l)  :  (2*j))  =  p'; 
vcombined ( 1 , j )  =  v; 
psidotcombined ( 1 , j )  =  psidot; 

%  If  attacker  is  already  neutralized  skip  to  next  one 
if  (NeutAtt ( j , DefenderTotal+1 )  ~=  0) 
continue; 

end 

for  k  =  1 : DefenderTotal 

%  Attacker  not  assigned  to  this  defender,  skip  to  next 
defender 

if  (  ( j  <Assign ( k, 1 ) )  |  |  ( j  >Assign ( k, 2 ) ) ) 
continue; 

end 
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%  Calculate  range  to  attacker 
Li  =  norm(p-Dpos(:,k)); 

Range ( k, i+1 , 1 )  =  t; 

Range ( k, i+1 , j +1 )  =  Li; 

Ai  =  (p-Dpos ( : , k) ) /Li ;  %  Unit  LOS  vector  to  attacker 

dVel  =  VelD ( : , k)  -  vel; 

Vc  =  Ai'*dVel;  %  Closing  velocity 

%  Calculate  LOS  angle 
temp  =  atan2 (Ai (2 ) , Ai ( 1 ) ) ; 

LOS  =  (temp- (pi/2-Dpsi (k) ) ) *180/pi; 
while  (abs(LOS)  >  180) 

LOS  =  -l*sign (LOS) * (360  -  abs(LOS)); 

end 

%  Sum  weighted  position  and  velocity  if  within  FOV 
if  (abs(LOS)  <=  Dfov/2) 
if  (Li  <  150) 

%  Attacker  within  FOV  and  attacking  range  of 
%  defender,  increase  counter  (one  time  step) 
NeutAtt(j,k)  =  NeutAtt(j,k)  +  1  ; 

%  If  attacker  stays  within  attacking  range  for  X 
%  secs  attacker  considered  neutralized 
if  (NeutAtt ( j , k) *dt  >  3) 

NeutAtt ( j , DefenderTotal+1 )  =  i; 
continue; 

end 

else 

NeutAtt (j,k)  =  0;  %  Reset  neutralize  counter  if 

attacker  goes  out  of  range 

end 

WPos ( : , k)  =  WPos ( : , k)  +  (Vc/Li)*p; 

WVel ( : , k)  =  WVel ( : , k)  +  (Vc/Li) *vel; 

W (k)  =  W (k)  +  (Vc/Li) ; 

else 

NeutAtt (j,k)  =  0;  %  Reset  neutralize  counter  if 

attacker  goes  out  of  FOV 

end 

end 

end 

%  If  all  attackers  neutralized,  end  simulation 
if  ( (prod (NeutAtt (:, DefenderTotal+1 ) )  ~=  0)  &&  (SimEnd  ==  0)) 

SimEnd  =  i; 

end 

%  Guidance  algorithm 
for  k  =  1 : DefenderTotal 
Defender ( i+1 ,  1)  =  t; 

Defender (i+1,  (2*k) : (2*k+l) )  =  Dpos ( : , k) ' ; 

%  Defender  lost  sight  of  attacker,  stop  guidance 
if  (abs(W(k))  <  le-15) 
continue; 

else  %  Attacker (s)  in  FOV 

%  Calculate  weighted  centroid  and  velocity 
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CPos  =  WPos ( : , k) /W (k) ; 

CVel  =  WVel ( : , k) /W (k) ; 

end 

Alog(i+l,  1)  =  t; 

Alog(i+l,  (2*k) : (2*k+l) )  =  CPos'; 

Li  =  norm (CPos  -  Dpos ( : , k) ) ;  %  Calculate  range  to  centroid 

Ai  =  (CPos-Dpos ( : , k) ) /Li ;  %  Unit  LOS  vector  to  centroid 

dVel  =  VelD(:,k)  -  CVel; 

Vc  =  Ai'*dVel;  %  Closing  velocity  of  centroid 

%  LOS  rate  of  centroid 

LOSrate  =  (dVel (1) *Ai (2) -dVel (2) *Ai (1) ) /Li; 
temp  =  atan2 (Ai (2 ) , Ai ( 1 ) ) ; 

LOS  =  (temp- (pi/2-Dpsi (k) ) ) ; 

%  Calculate  turn  rate  using  PN  guidance 
omega  =  6*L0Srate; 

%  Limit  turn  rate  by  MaxOmega 
if  (abs (omega)  >  MaxOmega) 

omega  =  sign (omega) *MaxOmega; 

end 

turn  =  -l*omega*dt;  %  Turn  angle 

Dpsi(k)  =  Dpsi(k)  +  turn;  %  Update  heading 

VelD ( : , k)  =  Dv* [sin (Dpsi (k) ) ; cos (Dpsi (k) ) ] ;  %  Update  velocity 

Dpos ( : , k)  =  Dpos ( : , k)  +  VelD ( : , k) *dt;  %  Update  position 

HeadingLog ( i+1 , 1 )  =  t; 

HeadingLog ( i  +  1 ,  k+1)  =  pi/2-Dpsi (k) ; 

end 

V  max  =  max(v,v  max); 

V  min  =  min(v,v  max); 

p  swarm(i+l,:)  =  [t  pcombined] ; 

V  swarm(i+l,:)  =  [t  vcombined] ; 
psi_dot ( i  +  1 ,  : )  =  [t  psidotcombined] ; 

%  Update  centroid  of  position  of  attackers 
APos  =  zeros (2, DefenderTotal) ; 
for  j =1 ; DefenderTotal 

for  k  =  Assign (j , 1)  rAssign (j , 2) 

APos ( : , j )  =  APos ( : , j )  +  pcombined ( 1 , 2 *k-l : 2 *k) '; 

end 

APos ( : , j )  =  APos (:, j )/ (Assign (j , 2) -Assign (j , 1) +1) ; 

end 

p_test (i+1, : )  =  p' ; 

end 
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g,_ _ _ _ _ _ _  _ _ _  _ _ _  _ _ _ _ _ _ 

o 

%  File:  PlotOneManyBadGuysV12 .m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  The  unit  LOS  vector  to  every  attacker  in  the  defenders'  FOV  is 
%  -  combined  to  obtain  the  heading  required  to  intercept  the  attackers. 
%  -  The  LOS  vectors  are  weighted  by  the  attackers'  range  and  velocity. 

%  -  This  version  divides  up  the  attackers  and  assigns  them  to  each 
%  -  defender.  The  defenders  will  stop  once  they  lose  sight  of  attackers 
%  -  assigned  to  them. 

%  -  This  version  also  includes  a  neutralizing  algorithm,  in  which  an 
%  -  attacker  is  neutralized  when  it  stays  within  the  FOV  and  attacking 
%  -  range  of  a  defender  for  a  set  amount  of  time. 

%  Inputs:  Initial  ship  and  swarm  parameters 
%  Outputs:  Defender  trajectory 

9-___  _  _____  ____ 

o 

global  vShip  psiShip  vO  vf  N  pO  DefenderTotal 
global  p  swarm  Defender  p  ship  Dindex  tf 

dt  =  tf/N;  %  Time  Step 
p_ship  =  zeros (N+1 , 3 ) ; 
for  i  =  0 : N 
t  =  i*dt; 

p_ship ( i+1 , : )  =  [t  pOship'  +  [vShip*cos (psiShip) 
vShip*sin (psiShip) ] *t] ; 
end 

%  Assigning  attackers  to  defenders 
%  This  algorithm  assumes  attackers  >=  defenders 

%  Attacker  groups  are  assigned  to  defenders  based  on  their  position 
Assign=zeros (DefenderTotal,  2)  ; 

Assign (DefenderTotal,  1)  =1; 

Assign (DefenderTotal, 2) =floor (BadGuyTotal/DefenderTotal) ; 
for  i=DefenderTotal-l : -1 : 1 

Assign ( i , 1 ) =Assign ( i+1 , 2 ) +1 ; 

Assign (i, 2) =Assign (i  +  1, 2) +floor (BadGuyTotal/DefenderTotal)  ; 

end 

for  i=l :mod (BadGuyTotal , DefenderTotal) 

Assign ( 1 : i , 2 ) =Assign ( 1 : i , 2 ) +1 ; 

end 

for  i=l :mod (BadGuyTotal, DefenderTotal) -1 
Assign ( 1 : i , 1 ) =Assign ( 1 : i , 1 ) +1 ; 

end 

%  Centroid  of  initial  position  of  attackers 
APos  =  zeros (2, DefenderTotal) ; 
for  j =1 : DefenderTotal 

for  i  =  Assign (j , 1) :Assign (j , 2) 

APos(:,j)  =  APos(:,j)  +  p0(:,i); 

end 

APos ( : , j )  =  APos (:, j )/ (Assign (j , 2) -Assign (j , 1) +1) ; 

end 
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%  Defenders  parameters 

MaxOmega  =1;  %  Max  turn  rate 

Dfov  =  120;  %  Field  of  view  (deg) 

Dv  =  25;  %  Speed  (m/s) 

%  Initialize  matrices 

range  max  =  0;  v  max  =  max(v0,vf);  v  min  =  min(v0,vf);  psidot  max  = 
p_test  =  zeros (N+1 , 2 ) ; 

p  swarm  =  zeros (N+1, (BadGuyTotal*2 ) +1 ) ;  v  swarm  = 
zeros (N+1 , BadGuyTotal+1 ) ;  psi_dot  =  zeros (N+1 , BadGuyTotal+1 ) ; 
pcombined  =  zeros ( 1 , BadGuyTotal*2 ) ;  vcombined  =  zeros ( 1 , BadGuyTotal ) 
psidotcombined  =  zeros ( 1 , BadGuyTotal ) ; 

Dpos  =  zeros (2, DefenderTotal) ;  Dpsi  =  zeros (DefenderTotal,  1)  ; 

Di  =  zeros (2  ,  DefenderTotal );  Defender  =  zeros (N+1,  DefenderTotal*2  +  l ) 
VelD  =  zeros (2, DefenderTotal) ; 

Range  =  zeros (DefenderTotal , N+1 , BadGuyTotal+1 ) ;  Dindex  = 
zeros (DefenderTotal,  1)  ; 

Alog  =  zeros (N+1,  DefenderTotal*2  +  l )  ;  HeadingLog  =  zeros (N+1, 
DefenderTotal+1 ) ; 

NeutAtt  =  zeros (BadGuyTotal ,  DefenderTotal  +  1); 

SimEnd  =  0; 

%  Initial  positions  and  heading  for  defenders 
for  i  =  1 : DefenderTotal 

Dpos(:,i)  =  pOship  +  (i-l)*[0;25]  +  [10;0]; 

%  Velocity  vector  pointed  towards  centroid  of  attackers 

VelD ( : , i )  =  Dv* (APos ( : , i ) -Dpos ( : , i ) ) /norm (APos ( : , i ) -Dpos ( : , i ) ) ; 

Dpsi(i)  =  pi/2  -  atan2 (VelD (2, i) ,  VelD(l,i)); 

end 

for  i  =  0 : N 
t  =  i*dt; 

Di  =  zeros (2,  DefenderTotal); 

DEN  =  zeros (DefenderTotal, 1) ; 
for  j  =  1 ; BadGuyTotal 

[p, V, psidot, vel]  =  MoreThanOneBadGuyData ( j  ,  tf ,  t)  ; 
pcombined (1,  ((2*j)-l)  :  (2*j))  =  p'; 
vcombined ( 1 , j )  =  v; 
psidotcombined ( 1 ,  j )  =  psidot; 

for  k  =  1 : DefenderTotal 

%  Attacker  not  assigned  to  this  defender,  skip  to  next 
%  defender 

if  (  ( j  <Assign ( k, 1 ) )  |  |  ( j  >Assign ( k, 2 ) ) ) 
continue; 

end 

%  Calculate  range  to  attacker 
Li  =  norm(p-Dpos(;,k)); 

Range ( k, i+1 , 1 )  =  t; 

Range ( k, i+1 , j +1 )  =  Li; 

Ai  =  (p-Dpos ( : , k) ) /Li ;  %  Unit  LOS  vector  to  attacker 

dVel  =  VelD ( : , k)  -  vel; 

Vc  =  Ai'*dVel;  %  Closing  velocity 
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%  Calculate  LOS  angle 
temp  =  atan2 (Ai (2 ) , Ai ( 1 ) ) ; 

LOS  =  (temp- (pi/2-Dpsi (k) ) ) *180/pi; 
while  (abs(LOS)  >  180) 

LOS  =  -l*sign (LOS) * (360  -  abs(LOS)); 

end 

%  Sum  weighted  LOS  vector  if  within  FOV 
if  (abs(LOS)  <=  Dfov/2) 
if  (Li  <  150) 

%  Attacker  within  FOV  and  attacking  range  of 
%  defender,  increase  counter  (one  time  step) 
NeutAtt(j,k)  =  NeutAtt(j,k)  +  1  ; 

%  If  attacker  stays  within  attacking  range  for  X 
%  secs  attacker  considered  neutralized 
if  (NeutAtt ( j , k) *dt  >  3) 

NeutAtt ( j , DefenderTotal+1 )  =  i; 
continue; 

end 

else 

NeutAtt (j,k)  =  0;  %  Reset  neutralize  counter  if 

attacker  goes  out  of  range 

end 

Di(:,k)  =  Di(:,k)  +  (Vc/Li)*Ai; 

DEN(k)  =  DEN(k)  +  (Vc/Li); 

end 

end 

end 

%  If  all  attackers  neutralized,  end  simulation 
if  ( (prod (NeutAtt (:, DefenderTotal+1 ) )  ~=  0)  &&  (SimEnd  ==  0)) 

SimEnd  =  i; 

end 

%  Guidance  algorithm 
for  k  =  1 : Def enderTotal 
Defender ( i+1 ,  1)  =  t; 

Defender (i+1,  (2*k) : (2*k+l) )  =  Dpos ( : , k) ' ; 

if  (abs(DEN(k))  <  le-15) 
continue; 

else  %  Attacker (s)  in  FOV 

Di(:,k)  =  Di(:,k)/DEN(k); 

Di  (  :  ,  k)  =  Di  (  :  ,  k) /norm  (Di  (  :  ,  k)  )  ; 

end 

turn  =  pi/2-atan2 (Di (2 , k) ,  Di(l,k))  -  Dpsi(k); 
if  (abs (turn)  >  MaxOmega*dt) 

turn  =  sign (turn) *MaxOmega*dt; 

Dpsi(k)  =  Dpsi(k)  +  turn; 

Di (  :  ,  k)  =  [sin (Dpsi (k) ) ; cos (Dpsi (k) ) ] ; 

else 

Dpsi(k)  =  pi/2-atan2 (Di (2 , k) ,  Di(l,k)); 

end 

Dpos ( : , k)  =  Dpos ( : , k)  +  Dv*Di ( : , k) *dt; 

HeadingLog ( i+1 , 1 )  =  t; 
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HeadingLog ( i+1 ,  k+l)  =  pi/2-Dpsi (k) ; 

end 

V  max  =  max(v,v  max); 

V  min  =  min(v,v  max); 

p  swarm(i+l,:)  =  [t  pcombined] ; 

V  swarm(i+l,:)  =  [t  vcombined] ; 
psi_dot ( i  +  1 ,  : )  =  [t  psidotcombined] ; 

p_test (i+1, : )  =  p' ; 

%  Update  centroid  of  position  of  attackers 
APos  =  zeros (2, DefenderTotal) ; 
for  j =1 : DefenderTotal 

for  k  =  Assign (j , 1)  rAssign (j , 2) 

APos ( : ,  j )  =  APos ( : , j )  +  pcombined ( 1 , 2 *k-l : 2 *k)  ' 

end 

APos ( : , j )  =  APos (:, j )/ (Assign (j , 2) -Assign (j , 1) +1) ; 

end 

end 
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C.  PERFORMANCE  RELATED  FUNCTIONS 


g, _ 

o 

%  File:  Probability_Function . m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  Calculate  the  probability  function  trajectories  x_a,x_d,  x_0 . 

%  -  Returns  a  cell  array  for  q  which  holds  the  probabilities  q{l}  for 
%  -  each  attacker's  survival  over  time.  Returns  an  array  p,  which  hold 
%  -  the  probability  at  each  time  point  for  HVU  survival.  Note  that  the 
%  -  final  objective  function  is  1-p. 

%  Inputs: 

%  Outputs: 

g, _  _ _ _ _ _ 

o 

function  [p,  q]  =  Probability  Function (x  a,  x  d,  x  0) 

global  CONSTANTS  OFFLINE_TRAJECTORIES  . . . 

PDF_VALUES  MESHED_PDF_VALUES . . . 

DISCRETIZATION_VALUES  MESHED_DISCRETIZATION_VALUES  . . . 
DIFFERENTIATION_MATRICES  . . . 

INTEGRATION_WEIGHTS  MESHED_INTEGRATION_WEIGHTS ; 

N  =  CONSTANTS. N; 

q  =  cell (1, CONSTANTS .ATTACKERS .Na) ; 
for  1  =  1 : CONSTANTS .ATTACKERS .Na 
q{ 1 }  =  zeros (N,  1 ) ; 
z  =  zeros (N,  1 ) ; 
for  i  =  1 :  N 
z  ( i )  =  0  ; 

attacker_positions  =  cell (1, CONSTANTS .ATTACKERS .Na) ; 
for  s  =  1 : CONSTANTS .ATTACKERS .Na 

attacker_positions { 1 , s }  =  x_a{s}(i,  1:2); 

%all  attacker  positions  influence  the  defender  hit  rate 

end 

point_to_evaluate  =  x_a{l}(i,  1:2); 
for  k  =  1 : CONSTANTS .Ns 

defender  position  =  x  d{k}(i,l:2); 
z(i)  =  z(i)+... 

feval (str2func (CONSTANTS . DEFENDER_HIT_RATE_FUNCTION) , 

k,  .  .  . 

point_to_evaluate, . . . 
defender  position, . . . 
defender  heading, . . . 
attacker_positions , . . . 
CONSTANTS .ATTACKERS .Na) ; 

end 

end 

for  i  =  1 : N 

%  This  line  of  code  assumes  that  the  integration  weights 
%  still  converge  over  the  partial  interval.  This  still 
%  needs  to  be  confirmed,  but  it's  definitely  true  for 
%  Euler's  method. 
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q{l}(i)  =  exp (-INTEGRATION_WEIGHTS{ 1 }  (1 :  i)  ' *z  (1 :  i) ) ; 

end 

end 

p  =  zeros (N,  1 ) ; 
z  =  zeros (N,  1 ) ; 
for  i=l : (N-1 ) 

for  1  =  1 : CONSTANTS .ATTACKERS .Na 

point_to_evaluate  =  x_0(i,l:2); 
attack;er_position  =  x_a{l}(i,  1:2); 

%only  the  position  of  the  particular  attacker  influences  the 
%attacker  hit  rate 
z(i)  =  z(i)+... 

q{l}  (i) *.  .  . 

feval (str2func (CONSTANTS .ATTACKER_HIT_RATE_FUNCTION) , . . . 

1,  .  .  . 

point_to_evaluate, . . . 
attacker_position) ; 

end 

end 

for  i  =  1 : N 

%  This  line  of  code  assumes  that  the  integration  weights 
%  still  converge  over  the  partial  interval.  This  still 
%  needs  to  be  confirmed,  but  it's  definitely  true  for 
%  Euler's  method. 

p(i)  =  exp (-INTEGRATION_WEIGHTS{ 1 }  (1 :  i)  '  *z  (1 :  i, 1) ) ; 

end 

g,  _  _  _  _  _  ___  ___  ___  ___  ___  ___  _ 

o 
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g,_ _ _ _ _ _ _ _ _ _  _ _ _  _ _ _  _  _ _ _ _ 

o 

%  File:  RunDingsSimulation . m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  Runs  Ding's  simulation  using  set  Heuristic,  Numberof Defenders ,  and 
%  -  Number  of  attackers.  pO,  the  array  of  starting  positions  of  the 
%  -  attackers,  is  also  assumed  to  be  passed  to  this  file  as  a  global 
%  -  variable  (it's  been  set  to  global  so  that  Kaminer  code  can  be  run 
%  -  without  editing) .  Returns  cell  arrays  x_a  and  x_d,  and  matrix  x_0, 
%  -  reformatted  to  be  in  the  form  used  in  the  optimal  control  set  up. 

%  -  Note:  Trajectories  still  need  to  be  interpolated  to  match  up  time 
%  -  points. 

%  Inputs:  Heuristic  function,  number  of  attackers,  number  of  defenders 
%  Outputs:  x_a,  x_d,  x_0 

g, _ _  _ _ _ _ _ _  _ _ _ _ _ _  _  _  _  _ _ 

o 

function  [x  a,  x  d,  x  0]  = 

RunDingsSimulation (Heuristic, NumberofDefenders, . . . 

Numberof Attackers ) 


global  N  DefenderTotal  BadGuyTotal  pO 
global  p  swarm  Defender  p  ship  Dindex  tf 

DefenderTotal  =  NumberofDefenders;  %  Total  number  of  defenders 
BadGuyTotal  =  Numberof Attackers ;  %  Total  number  of  attackers 

RunMoreThanOneBadGuy 

%  runs  some  of  the  Kaminer  code  with  whatever  value  pO  is  set  to. 

%  Note  that  this  version  of  the  Kaminer  code  has  a  preset  HVU 
%  trajectory  baked  into  it. 

feval (str2func (Heuristic)  )  ; 

%  runs  whichever  algorithm  is  set  as  Heuristic.  This  is  also  necessary 
%  for  running  the  rest  of  the  Kaminer  code. 

%  Trajectories  are  then  put  in  the  format  being  used  in  the  control 
%  implementation 

X  a=cell ( 1 , Numberof Attackers ) ;  %cell  array  of  attacker  trajectories 
for  j  =  1 : Numberof Attackers 
x_a { 1 ,  j } =zeros (N,  3); 

X  a { 1 , j } ( : , 3 ) =p  swarm ( 1 : N, 1 ) ;  %moves  time  to  the  third  column 
x_a { 1 , j }  ( : , 1 ) =p_s warm ( 1 : N , 2  *  j ) ; 
x_a { 1 , j } ( : , 2 ) =p_s warm (1:N,  2*j+l); 

end 

X  d=cell ( 1 , Numberof Defenders ) ;  %cell  array  of  defender  trajectories 
for  j  =  1 : Numberof Defenders 
x_d{ 1, j }=zeros (N,  3); 

X  d{ 1, j } ( : , 3) =Defender (1 :N, 1) ;  %moves  time  to  the  third  column 
x_d{ 1, j } ( : , 1) =De fender (1 :N, 2* j ) ; 
x_d{ 1, j } ( : , 2) =De fender (l:N,2*j+l) ; 

X  d{ 1, j } ( : , 4) =HeadingLog (1 :N, j+1) ;  %puts  heading  in  the  fourth 
column 

end 
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x_0=zeros (N,  3);  %HVU  trajectory 

X  0(:,3)=p  ship ( 1 : N,  1 ) ;  %moves  time  to  the  third  column 
x_0 ( : , 1 ) =p_ship ( 1 : N,  2 )  ; 
x_0 ( : , 2 ) =p_ship ( 1 : N,  3 )  ; 


9-___  __  _  _ 

o 

%  File:  Calibrate_Hit_Rate_Functions  . m 
%  Compiler:  MATLAB®  v7 . 1 0 . 0 . 4 9 9  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  Calibrates  hit  rate  functions  to  desired  shape 
%  Inputs:  None 
%  Outputs:  None 

g,  _  _  _  _  ___  _  _____  _______ 

o 

global  CONSTANTS 


% - Hit  Rate  Constants - % 

for  i=l :CONSTANTS .Ns 

CONSTANTS.DEFENDER_HIT_RATE{i} .alpha_theta  =  3; 

CONSTANTS. DEFENDER_HIT_RATE{ i } .beta_theta  =  3; 

%parameters  for  distribution  used  for  angle  effectiveness 
CONSTANTS. DEFENDER_HIT_RATE{ i } .mu_r  =  6; 
CONSTANTS.DEFENDER_HIT_RATE{i} .sigma_r  =  0.7; 

%parameters  for  distribution  for  radial  distance  effectiveness 
CONSTANTS. DEFENDER_HIT_RATE{ i } .max_angle  =  pi/3; 

%FOV  extends  plus  or  minus  this  angle  from  heading 
CONSTANTS.DEFENDER_HIT_RATE{i} .rho  =  1200; 

%radius  within  which  attackers  divide  defender  attention/rate 
CONSTANTS.DEFENDER_HIT_RATE{i} .sigma  =  .01; 

%tiny  number  for  standard  deviation  of  normal  cdf  that  smooth  rho 
x_star  =  exp (CONSTANTS. DEFENDER_HIT_RATE{ i } .mu_r  - 

(CONSTANTS. DEFENDER_HIT_RATE{i} .sigma_r) ^2) ; 
CONSTANTS. DEFENDER_HIT_RATE{i} . normalizing_constant_r  = 

1/lognpdf (x_star,  . . . 

CONSTANTS. DEFENDER_HIT_RATE{i} .mu_r, . . . 
CONSTANTS. DEFENDER_HIT_RATE{i} .sigma_r) ; 

end 

for  i=l :CONSTANTS .ATTACKERS .Na 

CONSTANTS.ATTACKER_HIT_RATE{i} .alpha  =  3; 

CONSTANTS. ATTACKER_HIT_RATE{i} .beta  =  18; 
CONSTANTS.ATTACKER_HIT_RATE{i} .cl  =  1; 

%changes  magnitude  of  rate  function 
CONSTANTS.ATTACKER_HIT_RATE{i} .c2  =  1/800; 

%inversely  changes  radius  of  rate  function 

end 

CONSTANTS . DEFENDER_HIT_RATE_FUNCTION  = 

' FOV_Def ender_Hit_Rate_logn_withdivider '  ; 
CONSTANTS .ATTACKER  HIT  RATE  FUNCTION  =  'Full  Beta  Attacker  Hit  Rate ' ; 
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g,_ _ _ _ _  _ _ _ _  _  _ _ _  _ _ _  _ _ _ 

o 

%  File:  Interpret_Results .m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  Inputs:  None 
%  Outputs:  None 

9-  __  ____________ 

o 

N_guess=500; 

Number o f Defender s=5; 

Numberof  Attack;ers=4  0 ; 

[x  a,  X  d,  X  0,  N  real ] =RunDingsSimulation (N  guess,  .  .  . 

Numberof Defenders , . . . 
Numberof Attackers ) 


colors  =  [ ' g ' ; ' r ' ; ' k ' ; ' c ' ; ' m ' ] ; 

plotters  =  [ ' xg ' ; ' xr ' ; ' xk ' ; ' xc ' ; ' xm ' ] ; 

figure;  hold  on 

for  j  =  1 : Numberof Attackers 

plot (x_a{ j }(:,!), x_a{ j }  ( : ,  2)  )  ; 

end 

plot (x_0 ( : , 1 ) , x_0 ( : , 2 ) ,  ' m ' ,  'LineWidth',  3 ,  ' DisplayName ' ,  ' HVU ' ) ; 

for  k  =  1 : Numberof Defenders 

plot (x_d{ k} ( : , 1 ) , x_d{ k} ( : , 2 ) ,  colors(k),  'LineWidth', 

3 ,' DisplayName ' ,  ['Defender  '  int2str (k) ] ) ; 
end 

hold  off 
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g,_ _ _ _ _ _ _ _  _  _  _ _ _  _ _ _  _ _ _ _ _ 

o 

%  File:  Plot_Hit_Rates . m 

%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Plot  hit  rate  functions 
%  Inputs:  None 
%  Outputs:  None 

g,  _  _  _  ____________ 

o 

global  CONSTANTS 

time  index  =  10; 

X  a  =  X  a  i ; 
x_d  =  x_d_a_i ; 
x_0  =  x_0_a_i ; 

figure 

colormap ('jet'); 
hold  on 

axis ( [-5, 8000, 0, 3200, 0, 1] ) 
view ( [-6  51] ) 
grid  on; 

spatial  incr  =  35;  %fineness  of  spatial  mesh 

[X,Y]  =  meshgrid ( -5 : spatial  incr : 8000, 0 : spatial  incr:3200); 

[i_length,  j_length] =size (X) ; 

for  s  =  1 : Numberof Attackers 
Z  =  0*X; 

for  i  =  l:i  length 

for  j  =  l:j  length 

point_to_evaluate  =  [X ( i ,  j )  ,  Y ( i ,  j ) ]  ; 
attacker  position  =  x  a{s}(time  index,  1:2); 

Z (i, j  )  = 

feval (str2func (CONSTANTS .ATTACKER_HIT_RATE_FUNCTION)  , 

s,  .  .  . 

point_to_evaluate, . . . 
attacker_position) ; 

end 

end 

mesh (X, Y, Z ) ; 

end 

for  s  =  1 : Numberof Defenders 
Z  =  0*X; 

for  i  =  l:i  length 

for  j  =  l:j  length 

point_to_evaluate  =  [X ( i , j ) , Y ( i , j ) ] ; 
defender  position  =  x  d{s}(time  index,  1:2); 
defender  heading  =  x  d{s}(time  index, 4); 
attacker^ositions  =~cell  (1,  CONSTANTS  .ATTACKERS  .Na)  ; 
for  a  =  1 : CONSTANTS .ATTACKERS .Na 

attacker_positions { a }  =  x_a{a} (time_index,  1:2); 

end 
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z (i, j  )  = 

feval (str2func (CONSTANTS . DEFENDER_HIT_RATE_FUNCTION) , 

s,  .  .  .  . 

point_to_evaluate, . . . 
defender  position, . . . 
defender  heading, . . . 
attack;er_positions ,  .  .  . 

CONSTANTS .ATTACKERS .Na) ; 

end 

end 

mesh (X, Y, Z ) ; 

end 


9-___  __  __  ___  ___  _  _  _  _ 

o 

for  s  =  1 : Numberof Attackers 

cline2D (x_a_i { s } ( : , 1 ) , x_a_i {s}(:,2),q{s}(:),' Cool ' ) ; 

end 

cline2D (x_0_a_i ( : , 1 ) , x_0_a_i ( : , 2 ) , p ( : ) , ' Autumn ' ) ; 

colormap ('jet'); 

for  s  =  1 ; Numberof Defenders 

plot (x_d_a_i { s }  ( : ,  1 )  ,  x_d_a_i { s }  ( : ,  2 )  ,  ' g '  ,  ' LineWidth ' , 3 ) ; 

end 

g,  _  _  _  _____  _  _  _  ___  _  _ 

o 
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g,_ _ _ _ _  _ _ _  _  _ _ _ 

o 

%  File:  Result_Plots .m 

%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 
%  64-bit  (win64) 

%  Description:  Plot  all  trajectories 
%  Inputs:  None 
%  Outputs:  None 

o  _  _  _____ 

o 


global  CONSTANTS  OFFLINE_TRAJECTORIES  . . . 

PDF_VALUES  MESHED_PDF_VALUES . . . 

DISCRETIZATION_VALUES  MESHED_DISCRETIZATION_VALUES  . . . 
DIFFERENTIATION_MATRICES  . . . 

INTEGRATION_WEIGHTS  MESHED_INTEGRATION_WEIGHTS ; 

NumberofAttackers  =  CONSTANTS .ATTACKERS .Na; 

Numberof Defenders  =  CONSTANTS . Ns ; 


—  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  — 

o 

%  Plot  a_i  trajectories 

g,  —  _  _  _  —  —  —  —  —  —  —  —  —  —  —  — 

o 

figure 
view (2 ) 
hold  on 

for  s  =  1 : Numberof Attackers 

cline2D (x_a_i { s }  ( : , 1 ) , x_a_i {s}(:,2),q{s}(:),' Cool '  )  ; 

end 

for  s  =  1 : Numberof Defenders 

plot (x_d_a_i { s }  ( : ,  1 )  ,  x_d_a_i { s }  ( : ,  2 )  ,  ' g ' ) ; 

end 

cline2D (x_0_a_i ( : , 1 ) , x_0_a_i ( : , 2 ) , p ( : ) , ' Autumn ' ) ; 

g, _ _ _ _ _ _ 

o 

figure 
hold  on 

plot (x_0_a_i ( : , 3 ) ,  p(:),'b') 
for  s  =  1 : Numberof Attackers 

plot (x_a_i { s }(:, 3) ,  q{  s  }(:),'  r  ')  ; 

end 
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%  File:  cline2D.m 

%  Compiler:  MATLAB®  v7 . 10 . 0 . 499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  This  function  plots  a  3D  line  (x,y,z)  encoded  with  scalar  color 
%  -  data  (c)  using  the  specified  colormap  (default=j et) ; 

%  -  SYNTAX:  h=cline2D (x, y, z, c, colormap) ; 

%  -  DBF  09/03/02 
%  Inputs: 

%  Outputs: 

g,  _  _  _  _  _  ___ 

o 

%  Edited  by  CLW,  2012,  to  make  2D,  and  to  make  spline  work  when  entries 
%  are  identical 

function  h=cline2D (x, y, c, cmap) ; 

if  nargin==0  %  Generate  sample  data... 
x=linspace (-10, 10,  101)  ; 
y=2  *x . ^2  +  3 ; 
z=sin(0.1*pi*x) ; 
c=exp ( z ) ; 
w=z-min ( z ) +1 ; 
cmap= 'jet'; 
elseif  nargin<3 

fprintf (' Insufficient  input  arguments \n ') ; 
return; 

elseif  nargin==3 
cmap= 'jet'; 
end 

cmap=colormap (cmap) ; %  Set  colormap 

%  Generate  range  of  color  indices  that  map  to  cmap 
yy=linspace (0,1, size (cmap, 1) ) ; 

cm  =  spline (yy, cmap ', c) ; %  Find  interpolated  colorvalue 
cm ( cm> 1 ) = 1 ; 
cm ( cm<  0 ) =  0 ; 

%  Lot  line  segment  with  appropriate  color  for  each  data  pair... 
for  i=l : length  (x) -1 

h (i) =line ( [x (i)  x (i  +  1) ]  ,  [y (i)  y(i  +  l)], 

'color', [cm(:,i)], ' LineWidth ' , 2 ) ; 

end 

return 
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g,_ _ _ _ _  _ _ _ _ _  _  _  _  _  _  _  _  _  _ _ _  _ 

o 

%  File:  Full_Beta_Attack;er_Hit  Rate.m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Defines  full  beta  attacker  hit  rate  function 
%  Inputs:  attacker_index,  point_to_evaluate,  attacker  position 
%  Outputs:  rate 

0,  _  _  _  ________  ___ 

o 

function  rate  =  Full  Beta  Attacker  Hit  Rate (attacker  index,  .  . . 

point_to_evaluate, . . . 
attacker_position) 


global  CONSTANTS 

alpha  =  CONSTANTS. ATTACKER_HIT_RATE{attacker_index} .alpha; 
beta  =  CONSTANTS .ATTACKER_HIT_RATE{attacker_index} .beta; 
cl  =  CONSTANTS. ATTACKER_HIT_RATE{attacker_index} .cl; 
c2  =  CONSTANTS. ATTACKER_HIT_RATE{attacker_index} .c2; 

X  star  =  (alpha-1)/ (alpha+beta-2 ) ;  %point  of  function  maximization 
normalizing  constant  =  l/betapdf(x  star,  alpha,  beta); 
distance  =  c2*norm (attacker  position-point  to  evaluate); 

if  distance  <  1 

rate  =  cl*normalizing  constant*betapdf (distance,  alpha,  beta); 

else 

rate  =  0; 

end 
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g,_ _ _ _ _ _ _ _  _  _ _ _  _  _ _ _ 

o 

%  File:  FOV_Defender_Hit_Rate  logn_withdivider 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Defines  defender  hit  rate  function 
%  Inputs:  defender  index,  point_to_evaluate,  defender  position, 

%  defender_heading,  attacker  positions,  Numberof Attackers 

%  Outputs:  rate 

q, _ _ _ _ 

o 

function  rate  =  FOV  Defender  Hit  Rate  logn  withdivider 

(defender  index, . . . 
point_to_evaluate, . . . 
defender  position, . . . 
defender  heading, . . . 
attacker_positions , . . . 
Numberof Attackers ) 


global  CONSTANTS 

alpha_theta  =  CONSTANTS . DEFENDER_HIT_RATE { defender_index } . alpha_theta; 
beta_theta  =  CONSTANTS . DEFENDER_HIT_RATE { defender_index } . beta_theta; 
mu_r  =  CONSTANTS. DEFENDER_HIT_RATE{defender_index} .mu_r; 
sigma  r  =  CONSTANTS . DEFENDER  HIT  RATE { defender  index}. sigma  r; 
max  angle  =  CONSTANTS . DEFENDER  HIT  RATE { defender  index}. max  angle; 
normalizing  constant  r  = 

CONSTANTS . DEFENDER  HIT  RATE { defender  index} .normalizing  constant  r; 
rho  =  CONSTANTS. DEFENDER_HIT_RATE { def ender_index } .rho;  ~  ~ 

sigma  =  CONSTANTS . DEFENDER_HIT_RATE { defender_index }. sigma; 

Attention  Dividing  Sum  =  1; 

%We  want  the  hit  rate  to  stay  constant  for  one  attacker,  but  decrease 
%proportionate  to  additional  nearby  attackers.  This  can  be  accomplished 
%by  the  following  shenanigans:  If  there  is  just  one  attacker  the 
%attention  dividing  sum  is  one.  if  more  attackers,  the  nearest  attacker 
%is  automatically  counted  as  the  1.  then  each  additional  attacker  is 
%added  on. 

if  Numberof Attackers>l 

%first  need  to  figure  out  which  attacker  is  nearest. 

distance  =  norm (attacker_positions {1,1} ( : ) -defender_position ( : ) ) ; 

index  of  min=l; 

for  i=2 : Numberof Attackers 

if  norm (attacker  positions { 1 , i }(:) - 

defender_position ( : ) ) <distance; 
distance  =  norm (attacker_positions { 1 ,  i }  ( : ) - 

defender_position ( : ) ) ; 

index  of  min  =  i; 

end 

end 

for  i=l : Numberof Attackers 
if  i~=index  of  min 

distance  =  norm (attacker_positions { 1 ,  i }  ( : ) - 

def ender_position ( : ) ) ; 

Attention  Dividing  Sum  =  Attention  Dividing  Sum  + 

normcdf (rho-distance,  0,  sigma); 
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end 


end 

end 

distance  =  norm (point  to  evaluate-defender  position); 

f  r  =  normalizing  constant  r*lognpdf (distance, mu  r,  sigma  r) ; 

angle  of  point  =  atan2 (point  to  evaluate (2 ) -defender  position (2 ),.. . 
point_to_evaluate ( 1 ) -defender_position ( 1 ) ) ; 

defender  heading  =  mod (defender  heading, 2 *pi ) ;  %mod  so  it's  easier 

%  adjust  heading  to  fit  into  atan2 ' s  format  for  polar  coordinates 
if  defender  heading>pi 

defender  heading=defender  heading-2*pi; 

end 

%angle_of_point 

9-5-9-5-9-9-9-5-9-5-9-5-9-9-Q-5-9-5-9-5-9-5-Q-5-9-9-9-9-9-9-9-5-9-5-9-5-9-5- 

oooooooooooooooooooooooooooooooooooooo 

%  There's  a  discontinuity  at  pi/-pi.  The  following  cases  deal 
%  with  that. 

oooooooooooooooooooooooooooooooooooooo 

if  pi-max  angle<=defender  heading  &&  defender  heading<=pi 

if  -pi<=angle  of  point  &&  angle  of  point<=-pi+max  angle 
angle  of  point  =  angle  of  point+2*pi; 

end 

end 

if  -pi<=defender  heading  &&  defender  heading<=-pi+max  angle 
if  pi-max  angle<=angle  of  point  &&  angle  of  point<=pi 
angle  of  point  =  angle  of  point-2*pi; 
end 
end 

9-5-Q-5-9-9-Q-5-9-5-Q-5-9-5-9-5-9-5-Q-5-9-5-Q-9-9-5-Q-5-9-9-Q-5-9-5-Q-5-9-5- 

oooooooooooooooooooooooooooooooooooooo 


delta_theta  =  (abs (defender_heading- 
angle  of  point) +max  angle)/ (2*max  angle); 
if  (delta_theta>0  &&  delta_theta<l ) 

x_star  =  (alpha_theta-l ) / (alpha_theta+beta_theta-2 ) ;  %point  of 
function  maximization 

normalizing  constant  =  l/betapdf(x  star,  alpha  theta,  beta  theta) 
f  theta  =  normalizing  constant*betapdf (delta  theta, alpha  theta, 
beta_theta) ; 
else 

f_theta=0 ; 

end 

rate  =  f  r*f  theta/Attention  Dividing  Sum; 
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o\°  o\°  o\° 


g,_  _  _  _  _  _  _  _ _ _  _ _  _  _ _  _  _ _ _  _ _ _  _  _  _  _  _ _ 

o 

%  File:  Calculate_Methods .m 
%  Compiler:  MATLAB®  v7 . 10 . 0 . 499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  This  file  is  problem  independent  and  should  never  have  to  be 
%  -  edited,  except  to  add  more  methods. 

%  -  This  function  creates  differentiation  weights  and  integration 
%  -  weights  for  a  given  method.  It  also  creates  the  arrays  for  the 
%  -  discretized  values  of  the  variables.  It  therefore  accesses  data 
%  -  sent  from  the  Problem  File,  by  using  the  global  CONSTANTS.  The 
%  -  weights  it  creates  are  problem  specific;  they  do  not  have  to  be 
%  -  renormalized  for  the  variable  intervals. 

%  Inputs:  Discretization,  Methods 
%  Outputs:  None 

g, _ _  _ _  _  _  _ _  _  ___  _ 

o 

%  ©  2012,  CLAIRE  WALTON.  Some  Rights  Reserved. 

0 _ 

O - 

%  Notes 

%  'Discretization'  should  be  an  array  with  the  values  for  how  fine 
%  the  discretization  is  for  time  and  parameter  space. 

%  E.g.  in  this  case  discretization  [5,9,9]  would  run  the  simulation  for 
%  5  time  steps  and  a  9x9  parameter  space. 

g, 

o 

%  'Methods'  should  be  an  array  with  the  same  dimension  as 
%  Discretization.  Each  entry  is  the  method  used  for  each  discretization 
%  variable. 

%  Method  0:  Euler 

%  Method  1 :  Pseudospectral  with  Igl  points 

0 _ 

O - 

function  Calculate  Methods (Discretization,  Methods) 

global  CONSTANTS  OFFLINE_TRAJECTORIES  . . . 

JOINT_PDF  MESHED_JOINT_PDF. . . 

DISCRETIZATION_VALUES  MESHED_DISCRETIZATION_VALUES  . . . 
DIFFERENTIATION_MATRICES  . . . 

INTEGRATION  WEIGHTS  MESHED  INTEGRATION  WEIGHTS 


Compute  details  for  Time  Domain 


0, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 

g, 

o 


g, 

o 

Note  % 

Currently,  the  time  domain  is  % 
being  calculated  with  N+1,  but  % 

parameter  space  is  using  N.  % 

This  was  inherited  from  previous  % 
code  and  will  probably  be  changed. % 

___  _  _  _  _  _  _  9- 

o 


if  Methods ( 1 ) ==0 

disp(' - ') 

disp ( ' Methods : ' ) ; 
disp ( ' Time :  Euler ' ) ; 
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[Dn,  Time  Array,  Weights] =euldiff (Discretization (1) ) ; 

end 

if  Methods ( 1 ) ==1 

disp(' - ') 

disp ( ' Methods : ' ) ; 

dispi'Time:  Pseudospectral  with  Legendre  Points'); 

[Dn,  Time  Array,  Weights] =legdiff (Discretization (1) ) ; 

%  Get  differentiation  matrix,  legendre  points,  and  quadrature  weights 
end 

% - Transformation - % 

%  Have  to  transform  these  from  the  % 

%  interval  [-1,1]  to  [T0,TF].  Equiv  % 

%  to  u  =  (1/2)  (TF-TO) t+ (1/2)  (TF+TO)  % 

g,  ___________  ____  9- 

o  o 

TO=CONSTANTS . Time . TO; 

TF=CONSTANTS . Time . TF; 

DISCRETIZATI0N_VALUES{1}  =  . 5 * ( TF+TO ) + . 5 * ( TF-T 0 ) . *Time_Array ; 

%  cell  array,  DISCRETIZATION_VALUES { 1 }  is  the  array  of  time  points 
DIFFERENTIATI0N_MATRICES{1}  =  ( 2 / ( TF-TO ) ) . *Dn; 

%  cell  array,  first  element  is  Dn  for  time. 

INTEGRATI0N_WEIGHTS{1}  =  . 5 * ( TF-TO ) . *Weights ; 

%  cell  array,  first  element  is  integration  weights  for  time 


%  Compute  details  for  Parameter  Space 


Size=CONSTANTS . ParameterSpace . Dimension+1 ; 

%  i.e.  length  (Discretization) 
for  i=2:Size 


if  Methods ( i ) ==0 

str= [' Parameter  ' , num2str (i-1) ,  Euler']; 
disp ( str ) ; 

[Dn,  Parameter  Array,  Weights] =euldiff (Discretization  (i) ) ; 

end 

if  Methods ( i ) ==1 

str= [' Parameter  ' , num2str (i-1)  , 

' :  Pseudospectral  with  Legendre  Points ' ] ; 

disp ( str ) ; 

[Dn,  Parameter  Array,  Weights] =legdiff (Discretization (i) )  ; 

%  Get  differentiation  matrix,  legendre  points,  and  quadrature  weights 
end 

WO=CONSTANTS. ParameterSpace. WO (i-1)  ; 

WF=CONSTANTS . ParameterSpace . WF ( i- 1 )  ; 

DISCRETIZATION_VALUES{i}  =  . 5 * ( WF+WO ) + . 5 * ( WF-WO ) . * Parameter_Array ; 
%  cell  array,  DISCRETIZATION_VALUES { 1 }  is  the  array  of  time  points 

DIFFERENTIATION_MATRICES{i}  =  ( 2 / ( WF-WO ) ) . *Dn; 

%  cell  array,  first  element  is  Dn  for  time. 

INTEGRATION_WEIGHTS{i}  =  . 5 * ( WF-WO ) . *Weights ; 

%  cell  array,  first  element  is  integration  weights  for  time 
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end 
disp  ( 


%%%%%Create  meshed  values  for  multiple  attack;ers%%%% 

%  there  is  no  mesh  for  the  time  domain,  but  it's  nice  to  keep  the 
%  dimension  numbers  matching  the  unmeshed  values 
MESHED_DISCRETIZATION_VALUES {!}=[]; 

MESHED_INTEGRATION_WEIGHTS {!}=[]; 

%  This  creates  the  meshes  for  evaluating  every  permutation  of 
%  discretization  values.  The  objective  function  still  needs  to  have  a 
%  nested  for  loop  iterating  through  the  two  meshes.  Each  i-th  column 
%  of  the  cell  is  the  value  for  the  i-th  attacker.  Each  row  is  a  unique 
%  permutation  of  Na  values.  Equivalent  meshes  are  created  for  the 
%  integration  weights  to  keep  track  of  which  weights  need  to  be  used 
%  for  each  permutation.  To  integrate,  take  the  product  over  all  Na 
%  columns. 

for  i=2:Size 

%  first  cell  arrays  are  used  to  trick  ndgrid  into  outputting  the  right 
%  dimension  grid 

discretization_values=cell (1, CONSTANTS .ATTACKERS .Na) ; 

integration_weights=cell (1, CONSTANTS .ATTACKERS .Na) ; 

if  CONSTANTS .ATTACKERS .Na>l 

[discretization_values { 1 , : } ]  = 

ndgrid (DISCRETIZATION_VALUES{i} ) ; 
[integration  weights { 1 ,:} ]  =  ndgrid ( INTEGRATION  WEIGHTS { i }) ; 

else 

discretization_values{l, : }  =  DISCRETIZATION_VALUES{i}; 
integration_weights{l, : }  =  INTEGRATION_WEIGHTS { i } ; 

end 

%  then  everything  is  reshaped  to  reduce  notational  confusion.  While 
%  the  meshgrid  outputs  Na  Na-dimensional  arrays,  these  can  be  indexed 
%  through  and  reread  as  Na  one-dimensional  vectors.  This  creates  a 
%  matrix  with  the  dimensions  (parameter_length''Na,  Na) 
MESHED_DISCRETIZATION_VALUES { i }  = 

zeros ( length (discretization_values { 1 ,  1 }  ( : ) ) ,  ... 

CONSTANTS .ATTACKERS .Na) ; 

MESHED_INTEGRATION_WEIGHTS{i}  = 

zeros (length (discretization_values { 1 ,  1 }  ( : ) ) ,  ... 

CONSTANTS .ATTACKERS .Na) ; 

for  s=l : CONSTANTS .ATTACKERS .Na 

MESHED_DISCRETIZATION_VALUES{i} (:,s)  = 

discretization_values { 1 , s } ( : ) ; 
MESHED_INTEGRATION_WEIGHTS{i} (:,s)  = 

integration_weights { 1 , s } ( : ) ; 

end 

end 
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g,_ _  ____  __  _ _  _ _ _  _ _ _ _ _ 

o 

%  File:  euldiff.m 

%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  This  function  calculates  the  differentiation  matrix  Dn  that  is 
%  -  obtained  by  using  forward  Euler's  method  on  n  equally  spaced 
%  -  points.  These  values  are  calculated  over  the  interval  [-1,1]  and 
%  -  need  to  be  transformed  to  the  interval  of  the  problem.  Note  that 
%  -  this  does  not  return  a  square  matrix.  It  returns  a  (N-l)xN  matrix 
%  Inputs:  N 

%  Outputs:  Dn,  x,  w 

0, 

o 

%  ©  2012,  CLAIRE  WALTON.  Some  Rights  Reserved. 

0 _ 

O - 

function  [Dn, x, w] =euldif f (N) ; 
h=2/ (N-1) ; 

X  =  (-l:h:l)  '  ; 

w  =  ones (N, 1 ) . *h; 

w (N, 1) =0; 

vl=ones ( 1 , N-1 ) ; 

v2=-l*ones (1,  N) ; 

temp= (diag (vl , 1 ) +diag (v2 ) ) . /h; 

Dn  =  temp ( 1 : N-1 , : ) ; 
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%  File:  legdiff.m 

%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  This  function  calculates  the  differentiation  matrix  Dn  that  is 
%  -  obtained  by  differentiating  the  Lagrange  Polynomials  at  the 
%  -  Legendre-Gauss-Lobatto  (LGL)  points.  It's  zero  on  the  main 
%  -  diagonal  except  at  l=k=l,  where  Dn(l,l)=  n(n+l)/4;  and  at  l=k=n; 
%  -  where  Dn (n, n) =-n (n+1 ) /4 .  n=  no  of  LGL  points.  For  the  other  LGL 
%  -  points  1  (~=)k,  we  have  Dn(l,k)=  Ln (xl) /Ln (xk) * (1/xl-xk) . 

%  -  This  routine  is  part  of  DIDO  Version  0.1 
%  -  Written  by  Fariba  Fahroo,  edited  by  I.  Michael  Ross 
%  -  Major  subfunctions  written  by  W.  Gragg 
%  -  Naval  Postgraduate  School,  Monterey,  CA  93943 
%  Inputs:  n 
%  Outputs:  Dn,  x,  w 

0,  _  _  _  ______  ___ 

o 

function  [Dn, x, w] =legdif f (n) ; 

[x  w] =lobatto (n) ; 
x=sort (x) ; 

%  initialize  Dn 
Dn=zeros (n)  ; 

Dn(l, 1)=- (n-1) *n/4; 

Dn (n, n) =n* (n-1 ) / 4  ; 

%  Calculate  the  legendre  polynomials  at  xi . 
p=0*eye (n) ; 

for  i=l:n;  s=x(i);  p(i,l)=l;  p(i,2)=s; 

for  j=2:n-l;  p ( i , j +1 )  =  ( (2  * j -1 ) * s*p ( i , j ) - ( j -1 ) *p  ( i , j -1 ) ) / j ;  end;  end; 

%  Fill  out  the  rest  of  matrix  Dn . 
for  1=1 :n;  for  k=l:n; 
if  l~=k, 

Dn(l,k)=p(l,n)/(p(k,n)*(x(l)-x(k))); 

end; 

end; end; 


0 _ 

O - 

function  [x,w]  =  lobatto (n, a, b) 

%  [x  w]  =  lobatto (n)  or  [x  w]  =  lobatto (n, alpha, beta) : 

q, 

o 

%  Computes  abscissa  and  weights  for  the  n-point  Gauss- Jacobi-Lobatto 
%  quadrature  formula  using  the  method  of  Gene  H.  Golub,  Some  modified 
%  matrix  eigenvalue  problems,  SIAM  Rev.  15  (1973)  318-334.  Another 
%  early  algorithm  for  this  is  by  David  Galant,  An  implementation  of 
%  Christoff el ' s  formula  in  the  theory  of  orthogonal  polynomials.  Math. 
%  Comp.  25  (1971)  111-113.  All  such  algorithms  should  be  "reviewed", 

%  in  light  of  recent  improvements  in  tqr  and  Cholesky  LR  algorithms. 

%  But,  this  algorithm  "ain't  bad". 

%  Copyright  (c)  23  August  1997  by  Bill  Gragg.  All  rights  reserved. 

%  lobatto  calls  mxt,  mxtj  and  tqr. 


73 


begin  lobatto 


if 

nargin  <  2 

a  =  0 ; 

b 

=  0; 

end 

m  = 

2'^  (a  +  b  + 

1) 

*beta (a+1 , b+1 ) ; 

us 

=  a  ==  b; 

n  = 

\ — 1 

1 

[a  b]  =  mxt  j  (n,  a,  b)  ; 

T  = 

mxt 

( a ,  b )  ; 

I  = 

eye (n) ; 

e  =  zeros (n, 1 ) ; 

e  (n) 

1; 

c  = 

(T  +  I) \e; 

c  =  c(n); 

d  = 

(T 

-  I)\e; 

d  = 

d  (n)  ; 

e  =  c  -  d ; 

c  = 

(c 

+  d) /e; 

d  = 

sqrt (2/e) ; 

a(n+l)  =  c; 

b(n) 

d; 

[X 

u]  =  tqr(a. 

b)  ; 

u  =  u  '  ; 

w  = 

m*u 

."2; 

%  "Purify"  formulas  in  the  ultraspherical  case, 
if  us 

X  =  (x  -  f  lipud  (x)  ) /2 ;  w  =  (w  +  f  lipud  (w)  ) /2 ; 

end 


%  end  lobatto 

q, _ _ _ 

o 

function  T  =  mxt(a,b,c) 

%  T  =  mxt(a,b,c)  or  T  =  mxt(a,b) : 

q, 

o 

%  T  =  mxt(a,b,c)  is  the  TRIDIAGONAL  MATRIX  with  diagonal  elements 
%  a(l:n),  subdiagonal  elements  b(l:n-l)  and  superdiagonal  elements 
%  c (1 :n-l)  . 

%  T  =  mxt(a,b)  is  the  HERMITIAN  tridiagonal  matrix  with  diagonal 
%  elements  real(a(l:n))  and  subdiagonal  elements  b(l:n-l). 

%  Copyright  (c)  1  December  1990  by  Bill  Gragg.  All  rights  reserved. 
%  Revised  21  November  1992. 


%  mxt  calls  no  extrinsic  functions. 


%  begin  mxt 
if  nargin  <  3 

a  =  real (a) ;  c  =  b ' ; 

end 


n  =  length (a) ; 
c  =  c ( 1 ; n-1 ) ; 


b  =  b (1 :n-l) ; 
z  =  zeros (n-1 , 1 ) ; 


if  n  <  500 

B  =  diag (b) ; 

C  =  diag (c) ; 

T  =  diag (a) ; 

else 

T  =  zeros (n) ; 


B  =  [  z  '  0 ;  B  z  ]  ; 

C  =  [z  C;  0  z '  ]  ; 
T  =  T  +  B  +  C; 


for  k 


1  :n-l 
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T(k,k;)  =  a(k;);  T(k;+l,k;)  =  b(k;);  T(k,k;+1)  =  c(k); 

end 

T(n,n)  =  a(n); 

end 

%  end  mxt 


g, _ _ _ _ _ _ _ _ _ 

o 

function  [a,b]  =  mxt j (n, alpha, beta) 

%  [a  b]  =  mxt  j  (n,  alpha,  beta)  ,  [a  b]  =  mxt  j  (n,  alpha)  ,  [a  b]  =  mxtj  (n)  , 

%  T  =  mxtj  (n,  alpha,  beta)  ,  T  =  mxtj  (n,  alpha)  or  T  =  mxtj  (n)  : 

q, 

o 

%  mxtj (n, alpha, beta) :  T  =  mxt(a,b)  is  the  Jacobi  matrix  whose 
%  characteristic  polynomial  p  is  (a  nonzero  scalar  multiple  of)  the  nth 
%  JACOBI  polynomial. 

%  The  eigenvalues  of  T  are  the  abscissas  of  the  nth  order  Gauss- 
%  Christoffel  quadrature  formula  for  the  weight  function  ((1  - 
%  t) ''alpha)  ((1  +  t) ''beta)  on  the  interval  -  1  <  t  <  1.  The  Gauss- 
%  Christoffel  weights  are  m(0)  times  the  squares  of  the  first  elements 
%  of  the  normalized  eigenvectors  of  T,  where  m(0)  =  b(0)''2  =  B  (alpha  + 

%  l,beta  +  1)2'' (alpha  +  beta  -  1)  is  the  total  mass. 

%  B  is  the  beta  function.  The  weight  function  is  positive  and 
%  integrable  if  alpha  +  1  >  0  and  beta  +  1  >  0. 

q, 

o 

%  mxtj (n, alpha)  takes  beta  =  alpha.  p  is  the  nth  ULTRASPHERICAL 
%  polynomial,  with  weight  function  (1  -  t''2) ''alpha  on  the  interval 
%  1  <  t  <  1 .  Special  cases  are  the  CHEBYSHEV  polynomial  of  the  FIRST 
%  KIND,  with  alpha  =  -  1/2,  and  of  the  SECOND  KIND,  with  alpha  =  1/2. 

q, 

o 

%  mxtj (n)  takes  alpha  =  beta  =0.  p  is  the  nth  LEGENDRE  polynomial, 

%  with  weight  function  w(t)  =  1  on  the  interval  -  1  <  t  <  1.  The 
%  quadrature  formula  here  is  originally  due  to  Gauss.  Christoffel 
%  generalized  Gauss'  formula  to  a  wide  class  of  weight  functions. 

%  Because  of  this  the  Gauss-Christof f el  weights  are  usually  called 
%  Christoffel  numbers. 

%  Copyright  (c)  2  February  1991  by  Bill  Gragg.  All  rights  reserved. 

%  mxtj  calls  mxt. 

%  begin  mxtj 
if  nargin  <  2 

alpha  =  0; 

end; 

if  nargin  <  3 

beta  =  alpha; 

end 

a  =  alpha;  b  =  beta;  c=a+b;  d=b-a; 

s(l)  =  d/(c  +  2);  t(l)  =  (a  +  l)*(b  +  l)/(c  +  2)''2/(c  +  3)  ; 


if  n  >  2 

d  =  c  *  d ; 
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o\o  o\o  fl)  o\o  o\o  o\o  o\o  o\o 


n=  (2:n)';  m=2*n;  mm  =  m  -  1 ;  mp  =  m  +  1 ; 

s (n)  =  d./(c  +  m)./(c  +  (m  -  2 ) ) ; 

t(n)  =  n.*(a  +  n).*(b  +  n).*(c  +  n)./(c  +  mm) . / 

(  (c  +  m)  . 


end 

a  =  s ( : )  ;  b  =  2*sqrt (t ( : ) ) ; 
if  nargout  <  2 

a  =  mxt (a, b) ; 

end 

%  end  mxtj 


^  2  )  .  / ( c  +  mp ) ; 


%  Problems . 

%  1.  Relate  T  =  mxt(a,b),  with  [a  b]  =  mxtj (n, 1/2),  with  the  negative 
%  second  difference  matrix  S  =  mxt(c,d),  with  [c  d]  =  mxs (n) . 

9~—  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  — 

o 

function  [lam,U]  =  tqr(a,b,U) 

%  [lam  u]  =  tqr(a,b)  or  [lam  U]  =  tqr(a,b,U) : 

g, 

o 

%  [lam  u]  =  tqr(a,b) : 

g, 

o 

%  The  column  lam  contains  the  eigenvalues  of  the  Hermitian  tridiagonal 
%  matrix  T  =  mxt(a,b)  computed  by  one  version  of  the  (real  symmetric) 

%  tqr  algorithm  with  Wilkinson's  shift.  The  column  u  contains  the 
%  first  elements  of  the  eigenvectors  of  T  normalized  to  be  nonnegative 

%  and  such  that  the  eigenvectors  are  unit  vectors.  In  practice  this  is 

%  an  0(n''2)  process.  If  u  is  omitted  only  the  eigenvalues  are 
%  computed.  The  computed  eigenvalues  are  real  and  are  sorted  to  be 
%  nonincreasing. 

g, 

o 

%  [lam  U]  =  tqr(a,b,U) : 

g, 

o 

%  This  replaces  the  input  U  by  UV  with  V  a  matrix  of  orthonormal  eigen- 

%  vectors  of  T.  If  the  input  U  is  I  the  output  U  is  V.  If  the  input  U 

%  is  unitary  with  AU  =  UT  then  the  output  U  is  unitary  with  AU  =  UD  and 
%  D  =  diag (lam) . 

g, 

o 

%  If  the  input  U  is  e(l) '  the  output  U  is  u' .  If  the  input  U  is 
%  [e(l) e (n) ']  the  output  U  is  [u';  v']  with  v  the  column  of  last 
%  elements  of  the  normalized  eigenvectors.  If  the  subdiagonal  elements 
%  of  T  are  all  nonzero  then  the  elements  of  v  alternate  in  sign,  at 
%  least  mathematically. 


Copyright  (c)  2  February  1991  by  Bill  Gragg.  All  rights  reserved. 
Revised  15  July  1994. 


tqr  calls  sgn. 
begin  tqr 


Ensure  that  T  is  Hermitian  and  shift  b  down  one  unit. 

=  real (a) ;  n  =  length (a);  b=  [0;  b(:)];  b=b(l:n); 

Initialize  U  if  required  and  execute  a  diagonal  unitary  similarity 
transformation  to  make  T  have  nonnegative  subdiagonal  elements. 

if  nargout  >  1 
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if  nargin  <  3 

U  =  zeros(l,n);  U(l)  =  1; 

end 

u  =  sgn(b);  u  =  cumprod (u) ;  U  =  U*diag(u); 

end 


b  =  abs (b) ; 

%  Scale  the  matrix  up  by  a  power  of  two  to  give  nearly  the  widest 
%  possible  exponent  range. 

scale  =  norm([a;  b*sqrt (2 ) ] ) ; 
scale  =  2''(1024  -  ceil  (log2  (scale)  ))  ; 
a  =  a*scale; 
b  =  b*scale; 

format  compact  %  Temporary  statements 

maxscale  =  max (abs ([a;  b] ) ) ;  %  for  display. 

%  "Do  tqr". 
for  m  =  n : -1 : 1 

%  Compute  the  mth  eigenvalue. 

for  its  =  0:10*n  %  its  is  the  iteration  index. 


%  Split  the  matrix  if  possible.  This  is  also  the  termination  test, 
for  k  =  m: -1 : 1 


if  k  >  1 

tol  =  abs(a(k-l))  +  abs(a(k)); 


end 

end 


if  tol  +  b(k)  ==  tol 

b(k)  =0;  break 

end 


if  k  ==  m 

break  %  b (m)  =  0.  a (m)  is  an  eigenvalue. 

else 


if  its  ==  10*n 

error ('tqr  iteration  did  not  terminate  in  lOn  steps 

end 


%  Compute  Wilkinson's  shift  w  as  a  perturbation  of  the 
%  Rayleigh  shift  r  =  a (m) .  As  the  algorithm  converges 
%  c  =  b (m)  -->  0. 

r  =  a (m) ;  c  =  b (m) ;  d  =  (r  -  a(m-l))/2;  s  =  abs (d) ; 
if  c  <  s 

s  =  c/s;  t  =  1  +  sqrt(l  +  s*s);  t  =  c*s/t;  %  t  <  c; 

else 
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s  =  s/c;  t  =  s  +  sqrt(l  +  s*s);  t  =  c/t;  %  t  <  c; 

end 

if  d  >  0 

w  =  r  +  t ; 

else 

w  =  r  -  t ; 

end 

%  Take  a  step  of  the  tqr  algorithm.  There  are  many  ways  to 
%  implement  the  inner  loop.  We  recently  found  the  fastest 
%  known  stable  form  in  terms  of  flops.  The  form  given  here 
%  is  elegant. 

c=l;  s=0;  p=w-a(k);  t=p; 

for  j  =  k:m-l 

%  Compute  the  two  by  two  reflector  stably  and  update  b(j) . 


oldc  =  c;  oldt  = 

t ;  q  =  b  ( j  + 1 )  ; 

u  =  abs (p) ; 

if  q  <  u 

V  =  q/u; 

r  =  sqrt(l  +  v*v) ; 

b(j)  =  u*r*s; 

u  =  sgn (p) ; 

c  =  u/r; 

s  =  v/r; 

else 

V  =  p/q; 

r  =  sqrt(l  +  v*v) ; 

b(j)  =  q*r*s; 

u  =  1; 

c  =  v/r; 

s  =  u/r; 

end 

%  Update  p,  t,  a(j) 

,  and  U(:,j:j+1)  if 

required . 

p  =  c* (w  -  a ( j+1) ) 

-  s*q*oldc;  t  =  c 

*p; 

a(j)  =  a(j+l)  +  t  - 

■  oldt; 

if  nargout  >  1 

i  =  j  :  j+1; 

U(:,i)  =  U(:,i)*[-c 

s ;  s  c  ]  ; 

end 

end 

%  Update  b (m)  ,  a (m) ,  and  U(:,m)  if  required, 
b (m)  =  abs(p)*s;  a (m)  =  w  -  t;  c  =  sgn(p); 

if  nargout  >  1 

U ( : ,m)  =  -  U ( : ,m) *c; 

end 

end 

end 

end 


%  Sort  and  prepare  the  output. 

[a  p]  =  sort (-a);  lam  =  -  a/scale; 


if  nargout  >  1 

U  =  U(:,p)  ; 


u  =  U  ( 1 ,  : )  '  ; 
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o\o  o\o 


if  nargin  <  3 

u  =  abs  (u)  ;  U  =  u '  ; 

else 

u  =  sgn(u);  U  =  U*diag(u'); 

end 

end 

%  end  tqr 


g,_ _ _ _ _ _ _ _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _ 

o 

function  W  =  sgn(Zl,Z2) 

%  W  =  sgn(Z)  or  W  =  sgn(Zl,Z2) : 

q, 

o 

%  For  z  a  complex  number  we  define  sgn  z,  the  SIGNUM  of  z,  as  z/|z|  if 
%  z  ~=  0  and  +  1  if  z  =  0.  Thus  sgn  z  is  the  same  as  matlab's  sign  z 
%  except  when  z  =  0.  We  always  have  I sgn  z|  =1,  apart  from  rounding 
%  errors. 

q, 

o 

%  With  the  first  call  W  is  the  Schur  (elementwise)  sgn  function  of  the 
%  complex  matrix  Z.  With  the  second  call  we  have  W  =  | Z1 | . *sgn (Z2) ; 

Copyright  (c)  19  January  1991  by  Bill  Gragg.  All  rights  reserved. 
Revised  29  May  1996. 

%  sgn  calls  no  extrinsic  functions. 

%  begin  sgn 
if  nargin  <  2 

W  =  sign(Zl);  W  =  W  +  (W  ==  0); 

else 

W  =  sign(Z2);  W  =  W  +  (W  ==  0);  W  =  abs(Zl).*W; 

end 

%  end  sgn 

%  Total  flops  (scalar  case,  see  csgn) :  TBC 
%  Problem. 

%  1.  Compare  this  function  experimentally  with  csgn.  Compare  with 
%  regard  to  both  execution  time  and  numerical  stability.  Is  matlab 
%  computing  sign  correctly? 

0 _ 

O - 


79 


g,_ _ _  _ _  _  _ _ _ _  _  _ _ _  _ _ _  _ _ _  _  _ _ _ _ 

o 

%  File:  Run_Optimization .m 
%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  Run  Optimization  functions  as  a  module  that  organizes  an  snopt  call 
%  -  for  both  searching  and  herding  problem 
%  -  Runs  on  general  form  agreed  on  for  objective  function 
%  Inputs:  Problem_Info 

%  Outputs:  X,  F,  xmul,  Fmul,  INFO,  run_time 

g, _ _ _ _ _ _  _ _ _  _ _  _  _  _  _  _ _ _  _ _  _  _ _ 

o 

%  ©  2012,  CLAIRE  WALTON.  Some  Rights  Reserved. 

0 _ 

O - 

%  Notes 

%  1.  This  file  should  be  run  after  all  offline  traj ectories/pdf s  are 
%  built 

%  2.  this  file  will  compute  sparsity  patterns  and  run  snopt 
%  3.  it  will  also  split  up  the  results  into  cost,  control,  and  searcher 
%  trajectories 

%  4.  run_time  will  measure  snopt  optimization,  and  will  not  include 

%  parsing  results 

g, 

o 

%  Claire  Walton,  01/19/12 

0 _ 

O - 

function  [x, F, xmul, Fmul, INFO,  run  time]  = 

Run  Optimization (Problem  Info) 

global  CONSTANTS  OFFLINE_TRAJECTORIES  JOINT_PDF 

DISCRETIZATION_VALUES  DIFFERENTIATION_MATRICES  INTEGRATION_WEIGHTS 

Examplel.spc  =  which (' Examplel . spc ') ; 
snspec  (Examplel.spc  ) ; 

[Objective  lower.  Objective  upper.  Dynamics  lower.  Dynamics  upper, .... 
Variables  lower.  Variables  upper]  = 
feval (str2func (Problem  Info . Optimization  Bounds)); 

xlow  =  Variables  lower; 
xupp  =  Variables_upper ; 

Flow  =  [Objective  lower;  Dynamics  lower]; 

Fupp  =  [Ob j ective_upper ;  Dynamics_upper ] ; 

%  F  is  the  vector  of  functions  made  by  concatenating  the  objective 
%  function  and  the  dynamics  constraints;  these  are  its  bounds. 

nF  =  length (Flow) ; 

%  dimension  of  the  function  vector  F 
n  =  length (xlow) ; 

%  number  of  state  variables 

X  =  feval (str2func (Problem  Info. Initial  Guess)); 

%x  =  zeros (n,l);  %this  is  for  no  initial  guess 
xstate  =  zeros (n,l); 

%  This  is  some  weird  thing  that  indicates  something  about  values  of  x. 

%  Setting  it  equal  to  zero  seems  to  be  related  to  no  initial  guess 
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Ob j Add  =  0; 

%  ObjAdd  is  a  constant  added  to  the  objective  function  for  printing 
%  purposes 

Ob j Row  =  1; 

%  ObjAdd  is  the  row  in  F  containing  the  objective  function.  We  will 
%  always  use  1 . 

xmul  =  zeros (n,l); 

%  This  has  something  to  do  with  the  vector  of  duals 

Fstate  =  zeros (nF, 1 ) ; 

Fmul  =  zeros (nF, 1 ) ; 

%  I  don't  know  what  these  are 

Start  =  1; 

%  The  value  of  Start  has  something  to  do  with  the  meaningfulness 
%  of  xstate  and  Fstate.  Start  =  1  knows  they're  not  meaningful. 

%  [iGfun, jGvar]  =  f eval ( str2func ( Problem_Inf o . Gradient_Sparsity) ) ; 

%  [iAfun, jAvar, A]  =  feval (str2func (Problem_Info . Linear_Gradient) ) ; 

disp(' - 

disp ( '  . 

disp('  .  running  snJac 

disp ( '  . 

[A, iAfun, jAvar , iGfun,  j Gvar]  = 
snJac ( Problem  Inf o . User Fun, x, xlow, xupp, nF) ; 

snseti  ('Derivative  option',  0); 
snseti  ('Major  Iteration  limit',  1000); 
snsetr  ('Major  feasibility  tolerance',  10'^  (-4)) 
snsetr  ('Minor  feasibility  tolerance',  10'^  (-5)  ) 
snsetr  ('Major  optimality  tolerance',  10^  (-4)) 
snsetr  ('Minor  optimality  tolerance',  lO^' (-5)  ) 

run  time  =  cputime; 

disp(' - 

disp ( '  . 

disp('  .  running  snopt 

disp ( '  . 

[x,  F,  xmul , Fmul , INFO]  =  snoptcmex (Start,  x,  xlow,  xupp,  xmul,  xstate. 

Flow,  Fupp,  Fmul,  Fstate,  ObjAdd,  ObjRow,  A,  .  .  . 
iAfun,  jAvar,  iGfun,  jGvar, . . . 

Problem  Inf o . UserFun) ; 

run  time  =  cputime  -  run  time; 

disp(' - '); 

disp ( '  Run  Time : ' ) ; 
disp (run_time) ; 

disp  (  ' - '  )  ; 


81 


g,_ _ _ _ _ _ _  _  _  _  _ _ _  _ _ _  _ _ _  _  _ _ _ _ 

o 

%  File:  Run_Problem.m 

%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description: 

%  -  This  file  is  problem  independent  and  should  never  have  to  be 
%  -  edited.  It  organizes  the  calls  to  three  other  functions: 

%  -  1.  The  function  to  create  the  differentiation  matrices  and 
%  integration  weights 

%  -  2.  The  function  to  build  any  attacker  or  hvu  trajectories  and  pdfs 
%  -  3.  The  function  the  organizes  the  snopt  call. 

%  -  The  purpose  of  the  file  is  mostly  to  create  the  global  variables. 

%  Inputs:  Problem_Info,  Discretization,  Methods 
%  Outputs:  Results,  build_time,  run_time,  INFO 

9-__  ___  ___  ___  ___  ___  ____________ 

o 

%  ©  2012,  CLAIRE  WALTON.  Some  Rights  Reserved. 

0 _ 

O - 

function  [Results,  build  time,  run  time,  INFO]  = 

Run  Problem ( Problem  Info, Discretization,  Methods) 

global  CONSTANTS  OFFLINE_TRAJECTORIES  . . . 

PDF_VALUES  MESHED_PDF_VALUES . . . 

DISCRETIZATION_VALUES  MESHED_DISCRETIZATION_VALUES  . . . 
DIFFERENTIATION_MATRICES  . . . 

INTEGRATION_WEIGHTS  MESHED_INTEGRATION_WEIGHTS 
build  time  =  cputime; 

H-i  '9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-9-S-'  ^  . 

'J.-LojJ  ^  OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO  )  f 

disp (Problem  Info. Name); 

disp(' - '); 

disp ( 'Discretization :  '  )  ; 
disp (Discretization)  ; 

Calculate_Methods (Discretization,  Methods) ; 

q, _ _ _ _ 

o 

%  Creates  DIFFERENTIATION_MATRICES ,  INTEGRATION_WEIGHTS ,  and 
%  DISCRETIZATION  VALUES. 

q, 

o 

%  DIFFERENTIATION_MATRICES  and  INTEGRATION_WEIGHTS  are  *cell  arrays* 

%  with  length  given  be  length  of  discretization  array.  Each  element  of 
%  the  cell  array  is  the  matrix/array  given  by  calculated  the 
%  differentiation  matrix/integration  array  of  the  corresponding 
%  discretization  variable. 

%  DISCRETIZATION_VALUES  is  also  a  cell  array. 

g, _ _ _ _ _ _  _ _ _ 

o 


feval (str2func (Problem  Info. Build  Simulation)); 

q, _ _ _ _ _ 

o 

%  Contributes  to  the  global  variable  TRAJECTORIES. 

%  TRAJECTORIES  may  include  TRAJECTORIES . HVU  and  TRAJECTORIES .ATTACKERS, 
%  depending  on  the  needs  of  the  objective  function  userfun. 

Q,  _  _  _  _  _  _  _  _  _  _  _  _  ___  _  ____ 

O 

feval (str2func (Problem  Info. Build  PDF)); 
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%  Creates  the  joint  PDF,  as  per  the  options  chosen  in  Problem  file. 
%  Put  results  in  global  cell  array  JOINT_PDF 

q, _ _ 

o 


build  time  =  cputime  -  build  time; 

disp(' - '); 

disp('  Build  Time:'); 
disp (build_time) ; 

disp(' - '); 

[x, F, xmul , Fmul , INFO,  run  time]  =  Run  Optimization (Problem  Info); 

g,  _  _  _  _  _  _  _  _  ___  __________________ 

o 

%  Organizes  the  snopt  call  and  calls  it.  Run  Optimization  will  access 
%  the  global  variables  created  here. 

q, _ _ _ _ _ _ _ _ _ _ _ 

o 


[Results]  =  ... 

feval (str2func (Problem  Inf o . Interpret  Results),  x, F, xmul , Fmul ) ; 


9-__  _  _  __  ___ 

o 

%  File:  Uniform  PDF 

%  Compiler:  MATLAB®  v7. 10. 0.499  (R2010a) 

%  64-bit  (win64) 

%  Description:  Forms  uniform  pdf 
%  Inputs:  discretization_values,  wO,  wf 
%  Outputs:  uniform  PDF 

q, _ _ _ 

o 

%  ©  2012,  CLAIRE  WALTON.  Some  Rights  Reserved. 

g _ 

o - 

function  PDF  =  Uniform  PDF (discretization  values , wO , wf) 
PDF  =  ones ( 1 , length (discretization  values )). /abs (wf-wO ) ; 
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